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ABSTRACT 


We  present  a  first-collision  model  for  the  evaluation  of  return  flux  from  the  exhaust  plume  of  a 
ring-symmetric  HF/DF  laser  in  LEO,  generated  by  an  incident  flux  of  ambient  molecules  traveling  at 
orbital  speed.  The  steady  plume  is  bounded  by  a  pair  of  lip-centered  rarefaction  fans,  and  unless 
spacecraft  attitude  enables  incident  air  molecules  to  reach  the  plume  through  the  cavitation  regions 
that  extend  beyond  these  fans,  the  spacecraft  is  shielded  from  ambient  scattering  by  its  own  plume. 
Assuming  hard-spheres  collisions,  the  first-collision  model  is  given  by  a  simple  closed-form  expression 
that  can  be  regarded  as  a  source  term  for  scattered  exhaust  molecules.  This  source  term  is  integrated 
numerically  throughout  the  fan,  yielding  the  flux  arriving  at  some  surface  "target  point". 
Quantitatively,  it  is  shown  that  for  a  typical  HF/DF  laser  exhaust  the  contamination  level  generated 
by  ambient  scattering  is  not  significant.  It  was  found  that  the  maximum  return  flux  of  HF  +  DF 
constitutes  about  2%  of  the  incident  ambient  flux;  this  ratio  will  be  nearly  constant  for  LEO 
altitudes.  The  value  of  this  flux  ratio  is  shown  to  be  dependent  on  the  molecular  collision  model;  it 
may  change  upon  replacing  the  hard-spheres  approximation  by  a  more  realistic  collision  model.  A 
possible  modification  of  spacecraft  charging  by  the  exhaust  was  examined,  including  production  of 
HF-  and  DF  — .  The  only  significant  effect  seemed  to  be  shadowing  of  the  downstream  half  of  the 
spacecraft  at  oblique  orbital  attitudes. 
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1.  INTRODUCTION 


This  presentation  is  part  of  a  study  on  the  gas  dynamics  of  ring- symmetric  exhaust  plumes  in 
space,  conducted  at  the  Naval  Postgraduate  School  in  Monterey.  A  ring-symmetric  jet  has  zero 
thrust,  which  makes  it  suitable  as  an  exhaust  configuration  for  various  open  loop  power  plants 
designed  to  produce  high  power  for  relatively  short  durations.  One  such  system  is  an  envisioned 
space-based  chemical  laser,  shown  schematically  in  Fig.  1-1.  In  the  case  of  a  chemical  laser,  a  ring- 
symmetric  configuration  would  also  enable  the  laser  radiation  to  emerge  in  the  form  of  an 
axisymmetric  beam. 

The  exhaust  nozzle  should  be  designed  to  bring  the  outgoing  flow  to  a  supersonic  speed  at  the 
nozzle  exit  surface.  The  near  field  of  a  free  jet  is  then  composed  of  an  inner  core  bounded  by  a  pair 
of  ring-symmetric  rarefaction  fans  centered  at  the  nozzle  lips  (Fig.  1-1).  Beyond  the  limiting 
characteristic  surface  of  the  centered  rarefaction  waves  (CRW),  a  near-vacuum  condition  prevails. 
For  the  purpose  of  continuum  gas  dynamic  analysis,  we  assume  it  is  a  perfect  vacuum. 

An  earth  orbiting  vehicle  is  subject  to  an  oncoming  stream  of  ambient  molecules  at  a  speed  of 
Ua*:8  (km/ sec),  in  a  direction  depending  upon  its  orientation  relative  to  the  orbital  velocity  vector. 
This  speed  is  sufficiently  high  to  cause  backscattering  of  exhaust  molecules  (see  schematic  description 
in  Fig.  1-2)  moving  at  speeds  appropriate  to  chemical  combustion  (about  2  to  4  km  s).  However, 
large  exhaust  plumes,  having  achieved  stationary  flow,  may  be  sufficiently  dense  at  their  outer  fringes 
to  effectively  trap  and  entrain  all  oncoming  ambient  molecules.  Thus,  ambient  scattering  may  be 
significant  only  in  selected  ranges  of  attitude  angles,  at  which  ambient  molecules  can  reach  the 
vicinity  of  the  spacecraft  by  traveling  almost  collisionlessK  through  cavitation  regions.  Exhaust 
molecules  that  may  be  "candidates"  for  ambient  scattering  will  hence  come  from  plume  segments 
flanked  by  cavitation  regions.  The  contribution  of  ambient  scattering  to  contamination  will  thus  be 
highly  dependent  upon  spacecraft  geometry  and  orientation.  This  may  well  affect  spacecraft  design 
and  operating  procedures. 

The  purpose  of  this  report  is  to  present  a  first-collision  model  for  estimating  the  flux  of  exhaust 
molecules  backscattered  from  the  fringes  of  the  plume  by  ambient  molecules,  along  with  results  of 
sample  flux  computations  performed  on  a  typical  HF  DF  laser  exhaust  configuration.  The  llow  field 
throughout  the  plume  is  assumed  to  be  governed  by  the  equations  of  continuum  gas  dynamics.  In 
principle,  the  »Lw  could  be  obtained  by  solving  the  governing  equations,  i.e..  the  equations  for 
stationary  isentropic  How  in  two-dimensional  axisymmetric  coordinates.  In  practice,  this  is  normally 


accomplished  by  integrating  the  flow  equations  in  characteristic  form,  using  some  finite  difference 
scheme  (method-of-characteristics).  We  have  pe.  formed  such  computations,  but  given  the  complexity 
of  applying  them  to  the  subsequent  integration  of  ambient  scattering  flux  (due  to  the  need  for  two- 
dimensional  interpolations  from  an  irregular  solution  grid),  we  opted  for  a  different  alternative  :  a 
closed-form  approximation  to  the  ring-symmetric  CRW,  based  on  an  analytic  expression  for  flow 
variables  along  characteristic  lines  that  fan  out  from  the  nozzle  lip. 

The  plan  of  this  report  is  as  follows.  In  Ch.  2  we  outline  the  approximation  to  the  ring- symmetric 
CRW  and  present  some  computation  results  that  demonstrate  its  accuracy.  In  Ch.  3  we  describe  the 
first-collision  model  and  the  3-D  spatial  integration  scheme  for  computing  the  flux  arriving  at  the 
cylindrical  spacecraft.  In  Ch.  4  some  results  of  backscattered  flux  of  corrosive  molecules  (HF+DF), 
showing  flux  variation  with  target  point  location  (Xs)  and  attitude  angles  (vj/A,  <pA)  are  presented.  In 
Ch.  5  we  take  up  the  subject  of  spacecraft  charging,  using  results  of  ambient  scattering  to  assess  the 
effect  of  laser  exhaust  on  spacecraft  charging.  This  is  followed  by  concluding  remarks  in  Ch.  6  and  a 
list  of  references  in  Ch.  7.  A  concise  description  of  the  flux  computation  code  "AMB"  is  given  in 
Appendix  A,  followed  by  the  code  listing. 
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2.  COMPUTATION  OF  THE  PLUME  FLOW  FIELD 
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Most  ambient  molecules  entering  the  CRW  that  flanks  the  exhaust  plume  are  stopped  within 
several  mean  free  paths  from  their  point  of  entry.  A  quantitative  estimate  of  ambient  back- scattering 
w'ould  thus  depend  on  the  flow  field  at  the  outer  (hypersonic)  fringes  of  the  lip-centered  CRW.  Even 
though  the  flow  in  those  regions  is  generally  past  the  surface  of  continuum  breakdown,  the  density 
there  is  reasonably  well  approximated  by  the  continuum  flow  field,  as  demonstrated  by  Bird's  Monte- 
Carlo  simulation  of  a  Prandtl-Meyer  expansion  to  vacuum  [1]  .  The  evaluation  of  ambient  scattering 
thus  calls  for  an  ancillary  computational  procedure  capable  of  rendering  the  continuum  flowr  field  at  a 
large  number  of  points  in  the  ring-symmetric  CRW  of  an  exhaust  plume.  This  method  uras  described 
in  a  recent  report  [2]  .  Here  wre  just  outline  the  key  ideas  and  main  results  of  this  approximation 
method. 

Our  analytic  approximation  to  a  ring-symmetric  CRW  is  formulated  as  follows.  In  a  planar  CRW 
(Prandtl-Meyer  flow')  all  flow  variables  are  uniform  along  the  characteristic  lines  that  fan  out  from  the 
corner  (we  assume  they  are  the  C+  family).  In  the  ring-symmetric  case  the  flow  near  the  corner 
approaches  asymptotically  a  corresponding  planar  CRW  flow,  which  we  term  the  associate  CRW. 
However,  the  gradients  along  C+  characteristics  at  the  corner  of  a  ring-symmetric  CRW  do  not 
vanish  as  in  a  planar  CRW.  The  key  idea  is  thus  :  evaluate  flow  gradients  in  C+  directions  at  the 
comer,  then  use  them  to  extrapolate  the  associate  CRW  along  C+  lines  to  a  finite  distance  from 
the  corner.  The  extrapolation  is  a  nonlinear  function  of  the  radial  coordinate  y  .  chosen  so  that  the 
ensuing  expression  conforms  exactly  to  the  flow  at  the  leading  (exit)  characteristic  C  +  (P,). 
Omitting  all  details  of  the  analysis,  the  resulting  approximation  is  presented  as  the  following  power- 


6(0. P) 


f(a.p)  =  f(0.p)  [y(«.p)/y(0.p)J 


where  f  is  the  streamtube  area  ratio  for  isentropic  flows  (  f=  1  at  a  sonic  point).  P  is  the  Mach 
number  of  a  particular  characteristic  line  at  the  corner,  a  is  a  coordinate  along  the  C  +  (P) 
characteristic  line  (  a  =  0  at  the  corner),  and  y  is  the  radial  coordinate  of  a  point  on  the 
characteristic  line  C  +  (P)  .  The  Mach  number  at  point  (a.p)  is  readily  determined  from  f(a.P)  using 
the  standard  relation  between  area  ratio  and  Mach  number  [3|  .  A  closed-form  expression  for  6(0. P) 
was  developed  but  is  not  given  here;  instead,  this  function  is  shown  in  Fig.  2-1.  We  note  that  6 
approaches  the  asymptotic  value  of  2/(3-y)  as  p  increases  to  infinity,  and  that  generally 
1  <  6(0. P)  <  2  so  that  streamtubes  diverge  at  a  rate  intermediate  between  that  of  cylindrical  and 
spherical  expansion  flows. 


/*  /•  «"•  V"  v  >  V%.'  V  V  V  */\*  V 


Clearly,  in  an  isentropic  flow  all  thermodynamic  variables,  and  in  particular  density,  can  be 
evaluated  from  f .  This  approximation  is  readily  applied  to  the  hypersonic  portions  of  a  ring- 
symmetric  CRW  since  it  turns  out  that  characteristic  lines  are  nearly  straight  there,  which  means  that 
the  characteristic  line  C  +  (p)  passing  through  a  given  point  can  be  readily  determined.  As  a 
demonstration  of  the  degree  of  accuracy  obtainable  from  this  approximation,  we  show  in  Fig.  2-2  the 
variation  of  Mach  number  along  a  characteristic  line  in  the  ring- symmetric  CRW,  compared  with  an 
accurate  method-of-characteristics  computation.  This  comparison  demonstrates  that  the  analytic 
approximation  is  reasonably  accurate  to  nearly  ten  comer-radii  away  from  the  corner. 
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3.  AMBIENT  SCATTERING 


When  a  rocket  or  laser  exhaust  is  released  into  space  from  an  earth-orbiting  spacecraft,  it 
encounters  an  oncoming  stream  of  ambient  molecules  flowing  at  the  orbital  speed  of  UA^8  (km  sec). 
At  altitudes  higher  than  200  (km),  the  air  air  mean  free  path  exceeds  250  (m),  so  that  it  is 
considerably  larger  than  almost  any  spacecraft.  Consequently,  ambient  molecules  would  hardly  be 
subjected  to  a  shock  transition  prior  to  their  impact  at  the  spacecraft  or  exhaust  plume.  In  this 
chapter  we  describe  the  formulation  of  the  first-collision  model  in  Section  3.1  and  then  proceed  to 
present  the  derivation  of  the  flux  integration  scheme  for  hard-sphere  collisions  in  Section  3.2. 

3.1  First  Collision  Model 

The  highest  ambient  number  density  that  we  consider  for  earth-orbiting  spacecrafts  is 
Na=  1  *  1016  (m*3),  which  roughly  corresponds  to  Sunspot  Maximum  at  200  (km)  [4]  .  The  typical 
laser  exhaust  (Table  4-i)  would  reach  a  number  density  of  about  2  x  1019  (m*3)  at  the  very  high  Mach 
number  of  30.  Hence,  ambient  flux  constitutes  just  a  slight  perturbation  to  the  near-field  portion  of 
a  typical  laser  exhaust  plume.  Obviously,  ambient  molecules  that  penetrate  the  plume,  would 
subsequently  be  entrained  by  the  main  flow.  But  how  far  do  they  penetrate?  And  would  exhaust 
molecules  scattered  by  them  reach  the  spacecraft?  In  seeking  answers  to  these  questions,  we  are  led 
to  some  interesting  observations  concerning  ambient  scattering. 

Consider  the  HF  laser  depicted  in  Fig.  1-1.  The  spacecraft  diameter  is  5  (m)  and  the  centrally 
located  ring- symmetric  nozzle  is  1  (m)  wide.  Typical  operating  conditions  (Table  4-1)  are  assumed. 
They  are  based  on  some  experimental  HF  DF  laser  studies  conducted  at  TRW  [5.6]  .  Suppose  that 
the  spacecraft  axis  is  normal  to  the  orbital  velocity  vector  (normal  incidence).  Let  the  plane  of 
incidence  be  the  plane  defined  by  the  intersection  of  the  spacecraft  axis  with  the  orbital  velocity 
vector.  The  probability  that  an  ambient  molecule  traveling  in  the  plane  of  incidence  would  reach  the 
spacecraft  collisionlessly  is  e.xp(  -  q),  where  q  is  its  expected  number  of  collisions  with  exhaust 
molecules.  We  define  the  number  q  as  "molecular  thickness",  in  analogy  to  "optical  thickness".  So  in 
order  to  determine  the  extent  to  which  ambient  molecules  at  normal  incidence  reach  the  spacecraft, 
we  seek  the  distribution  of  radial  molecular  thickness  as  function  of  distance  from  the  spacecraft  mid- 
plane  (normal  to  axis  at  its  midpoint). 


For  this  purpose  we  computed  the  ring-symmetric  exhaust  flow  field,  using  a  semi-inverse 
marching  characteristics  scheme  [7]  .  The  marching  was  in  the  radial  direction,  starting  with  uniform 
flow  at  the  nozzle  exit;  the  computation  was  carried  on  until  it  became  evident  that  even  at  a 
distance  of  20  (m)  from  the  mid-plane,  the  radial  molecular  thickness  was  well  over  40.  The  entire 
spacecraft  was  thus  shielded  from  any  ambient  scattering  at  (or  near)  normal  incidence.  This 
shielding  effect  has  two  significant  implications  which  we  discuss  briefly  below. 

(a)  It  is  present  only  during  stationary  exhaust  flow.  At  startup  and  shutdown  phases,  ambient 
scattering  may  be  substantial  even  at  normal  incidence. 

(b)  During  the  stationary  phase,  ambient  scattering  is  substantial  only  at  attitude  angles  that  enable 
ambient  molecules  to  reach  the  vicinity  of  the  plume  by  traveling  through  "molecularly  thin" 
cavitation  regions  that  flank  the  plume.  We  thus  anticipate  a  decisive  dependence  of  ambient 
scattering  on  attitude  variations,  whenever  those  variations  steer  the  spacecraft  into  or  out  of  a 
shielded  posture. 

As  a  first  attempt  at  a  quantitative  estimate  of  ambient  scattering  flux,  we  have  formulated  a 
simple  first-collision  model  of  this  effect.  In  the  sequel  we  present  an  outline  of  the  model,  along  with 
some  sampie  results  evaluated  for  an  HF  laser  configuration  identical  to  that  considered  for  the 
shielding  elfect  mentioned  above. 

The  basic  idea  is  the  following.  Ambient  molecules  entering  an  exhaust  plume,  require  several 
collisions  to  become  fully  "accommodated"  with  the  main  flow  (i.e..  to  be  entrained  by  the  main  flow 
at  the  prevailing  flow  velocity  and  temperature).  One  may  reasonably  approximate  this  process  by 
considering  just  one  collision  -  the  first. 

With  the  help  of  some  additional  assumptions,  we  were  able  to  derive  a  closed  form  expression  for 
the  flux  of  exhaust  molecules  that  arrive  at  the  spacecraft  following  a  first  collision  with  an  ambient 
molecule.  The  main  assumptions  of  this  model  are  : 

U)  FIRST  COLLISIONS:  Only  first  collisions  for  either  ambient  or  exhaust  molecules  are 
considered.  Hard-spheres  elastic  collisions  are  assumed.  Upon  a  second  collision  of  either  an 
ambient  or  an  exhaust  molecule,  it  is  considered  "lost"  (i.e..  it  joins  the  main  flow).  Collisions 
of  ambient  molecules  with  spacecraft  surfaces  are  ignored.  Ambient  molecules  are  assumed  to 
traverse  cavitation  regions  collisionlessly. 


(2)  COLD  FLOW:  The  oncoming  ambient  air  flow  is  deemed  "cold";  i.e.,  all  molecules  move  at 
the  uniform  orbital  velocity.  The  same  "cold"  assumption  is  applied  to  the  exhaust  flow,  since 
most  ambient  scattering  takes  place  at  plume  regions  of  very  high  Mach  numbers  (well  over  10, 
in  the  present  case). 

(3)  CRW  Flow  Field:  ring-symmetric  CRW  flow  field  is  determined  from  the  power-law 
approximation  described  in  Ch.  2  above.  This  approximation  approaches  Prandtl-Meyer  flow 
at  points  whose  distance  from  the  nozzle  lip  is  much  smaller  than  the  spacecraft  radius. 

Based  on  these  assumptions,  ambient  scattering  is  represented  as  a  source  term  for  side-scattered 
exhaust  molecules,  distributed  throughout  the  lip-centered  rarefaction  fan.  The  total  flux  arriving  at  a 
specified  point  on  the  cylindrical  spacecraft  is  readily  computed  by  integrating  numerically  that  source 
distribution  over  the  entire  ring-fan. 

The  highlights  of  the  spatial  integration  scheme  (Fig.  3-1)  are  as  follows.  The  limiting 
characteristic  surface  (M=  oo)  of  the  ring- symmetric  CRW  is  divided  into  surface  elements  formed  by 
dividing  the  surface  into  a  set  of  ring-strips  which  are  subdivided  in  the  circumferential  (azimuthal) 
direction  (<p)  into  surface  elements.  The  line-of-sight  (Cl)  from  the  "target  point"  on  the  spacecraft  to 
the  center  of  each  surface  element  is  extended  into  the  ring- symmetric  CRW,  and  flux  integration 
using  the  first-collision  source  term  with  appropriate  weight  factors  is  performed  along  this  line  until 
convergence  is  attained.  Contributions  from  each  surface  element  are  summed,  taking  care  to 
disregard  portions  of  the  ring-symmetric  CRW  that  are  shadowed  by  the  cylindrical  spacecraft  (either 
the  line-of-sight  or  the  trajectory  of  oncoming  ambient  molecules  may  be  shadowed).  Some  further 
details  of  the  flux  integration  scheme  and  hard-spheres  collisions  are  provided  in  Section  3.2  below. 


3.2  Flux  Integration  Scheme 

The  description  of  the  first  collision  model  is  hereby  supplemented  with  an  outline  of  the 
expressions  used  in  the  flux  integration  and  their  derivation.  The  integration  scheme  for  flux  arriving 
at  point  Xs  on  the  spacecraft  is  depicted  in  Fig.  3-1.  Note  that  only  the  plane  of  incidence  is  shown 
in  Fig.  3-1;  at  other  azimuth  angles  the  geometry  is  not  co-planar,  so  3-D  geometrical  expressions  arc 
used  ro  get  the  coordinates  (\j/.<p  and  radial  distancete  (y2  +  z2)^2)  from  12  and  S;  the  derivation  of 
these  geometrical  relations  is  straightforward,  so  that  we  omit  these  details  in  the  present  report.  The 
total  number  flux  Q;( Xs)  of  i  exhaust  molecules  arriving  at  point  Xs  is  given  by  the  following 
expression  : 
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Qj(Xs)  -  S  d3ft  cosas  L  J  dS  <rik  h;  N(S)  hk  NA  |  U(S)  -  UA  |  expl  -  nk(S)l  Pik(S,  -  ft)  expl  -  nik(S)l 
tcs) 

nk(S)  =  ^  |  dr  <rik  h;  N(to  |u(to-uA|/|uA|  (3-i) 

5  ^ 

nik(S)  =  L  J  dS'  ff;.  h.  N(SO  I  uik(S)  -  U(S0 1  / 1  uik(S)  I 
J  0 

( )j  ( )j  —  Exhaust  species  ( )k  -  Ambient  species 


These  expressions  are  interpreted  as  follows.  The  collision  depicted  in  Fig.  3-1  is  between  exhaust 
molecule  nij  and  ambient  molecule  mk.  The  exhaust  molar  fractions  hj  and  ambient  molar  fractions 
hR  are  assumed  uniformly  constant,  and  so  are  the  ambient  velocity  UA  and  number  density  NA.  The 
exhaust  velocity  U(S)  and  number  density  N(S)  are  function  of  the  location  in  the  flow  field  defined 
by  3  and  S.  These  flow  variables  are  computed  by  first  evaluating  the  coordinates  of  point  f£.S 
(Fig.  3-1)  in  the  ring-symmetric  CRW  from  the  3-D  geometry,  and  then  employing  the  power-law 
approximation  outlined  in  Ch.  2  above,  to  get  all  flow  variables  for  a  nng-symmetric  CRW.  In  this 
computation  we  exploit  the  fact  that  characteristic  lines  fanning  out  from  the  nozzle  lip  are  nearly 
straight  lines  at  the  low  pressure  side  of  the  ring- symmetric  CRW. 

The  ft  integration  is  performed  numerically  according  to  the  scheme  outlined  in  Section  3.1  above, 
as  a  summation  over  elements  of  solid  angle  ( A3ft )  subtended  by  area  elements  on  the  limiting 
characteristic  cone  (y  =  yf). 

The  S  integration  is  considerably  more  complex.  The  integrand  for  this  integration  is  derived  as 
follows.  Denote  by  L  the  line-of-sight  distance  between  point  Xs  and  fan  point  ft.S.  A  volume 
element  at  the  fan  point  is  given  by  Av  =  L'  AS  A^ft.  The  number  of  ik  pair  collisions  in  Av  per  unit 
time  is  <rjk  hj  N(S)  hk  N  v  |U(S)  -  U  J  e.\p[-  qk(S)J  Av  .  where  \(S)  denotes  the  expected  number  of 
collisions  of  ambient  molecule  k  with  any  exhaust  molecule,  between  its  point  of  entry  into  the  plume 
and  point  ft.S.  We  now  multiply  this  term  by  exp[  —  r|jk(S )|  which  is  the  probability  that  exhaust 
molecule  i  scattered  by  ambient  molecule  k  would  travel  from  point  ft.S  to  point  X  collisionlessly, 
where  Hjk(S)  is  the  expected  number  of  collisions  for  this  path  segment.  (Note  that  in  Eq.  (3-1)  the 
summation  in  the  expression  for  n,k(S)  is  over  all  exhaust  species  j). 


The  final  step  in  constructing  the  integrand  for  the  S  integration  involves  the  post-collision 
directional  distribution  function  pit(s,-n),  whose  derivation  will  be  given  in  the  sequel.  We 
multiply  the  integrand  by  P(k(S,  ~  it)  &3Clc  which  is  the  fraction  of  i  exhaust  molecules  scattered  by  k 
ambient  molecules  into  a  solid  angle  element  A3f2c  about  the  unit  vector  —  Cl.  Considering  the  flux 
arriving  at  a  surface  area  element  AAs  around  point  Xs,  the  solid  angle  element  subtended  by  AAs  is 
A 3Clc  =  AAs  cosas  /L2.  Eq.  (3-1)  for  Qj(Xs)  now  follows  upon  dividing  the  resulting  expression  by 
AAs,  thus  referring  the  arriving  flux  to  a  unit  area  at  the  point  of  arrival  Xs. 

Numerically,  the  S  integration  was  performed  using  the  classical  Runge-Kutta  scheme  (fourth 
order).  The  integration  for  Hjk(S)  and  Hk(S)  has  to  be  repeated  at  each  point  S.  We  found 
reasonable  convergence  with  4  points  in  the  Hk(S)  integration  and  6  points  in  the  azimuth  integration. 
The  S  integration  was  terminated  when  convergence  was  attained  (this  is  the  meaning  of  the  upper 
limit  20  in  the  S  integral  in  Eq.  (3-1)).  The  summation  over  new  strips  on  the  limiting  cone  (i|/  =  vj/f) 
was  also  terminated  upon  convergence.  The  CPU  time  consumed  per  target  point  was  about  100 
(sec)  on  IBM  3033  mainframe. 

We  now  take  up  the  derivation  of  an  expression  for  the  post-collision  directional  distribution 
function  p;k(S,  fi),  which  we  denote  hereafter  as  P(  ~  Cl).  We  adopt  the  pair-collision  notation 
presented  in  Fig.  3-2  for  the  hard-sphere  collision  analysis. 

As  a  consequence  of  conservation  of  momentum  and  energy  (elastic  collisions),  the  center-of-mass 
velocitv  C  and  the  magnitude  of  the  relative  velocity  C  are  unchanced  bv  the  collision  181  .  The 
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post-collision  velocities  are  given  by  : 


CJ  -  cm  +  n2c; 


C2  =  Cm  " 


Cr  “  Ct  -  C2 


c*  =  c*  -  c* 

r  '“l  v"> 


B,  =  m,/(m,  +  m-,) 


B-,  =  m^lm,  +  m,l 


Cm  -  l*|C|  +  fiCj 


C'  -  Cr 


The  only  free  parameter  in  the  expressions  for  post-collision  velocities  is  the  orientation  of  the 
post-collision  relative  velocity  C  .  This  orientation  is  uniformly  likely  to  be  in  any  direction  in  <pace 
when  hard-spheres  collision  is  assumed  [S]  ,  as  represented  by  the  spherical  scattering  envelope  in 


Fig.  3-3.  The  probability  of  obtaining  ^  in  solid  angle  element  A3H  about  -Q.  (Fig.  3-3)  is  given 
by  : 

P(-3)  =  (l/4«|p2Cr|2)(AA/A3n)  -  (l/47t|cos6|)(|c*|2/U2Cr|2)  (3-3) 

where  AA  is  an  area  element  on  the  scattering  envelope,  whose  projection  on  a  plane  normal  to  is 
AA  |cosS|.  We  note  that  the  origin  of  Cm  in  Fig.  3-3  is  external  to  the  scattering  envelope,  resulting 
in  two  possible  scattering  elements  on  the  sphere.  In  all  the  cases  that  we  computed,  however  (see 
Ch.  4  below),  that  point  was  found  to  be  always  internal,  so  that  there  was  only  a  single  scattering 
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solution  with  post-collision  velocity  Uik(S)  pointing  at  the  spacecraft  for  any  ik  pair  collision. 


4.  RESULTS  AND  DISCUSSION 


We  performed  several  computations  of  return  flux  generated  by  ambient  scattering,  aimed  at 
demonstrating  the  expected  flux  level  and  its  variation  with  spacecraft  target  point  and  orbital 
attitude  angles.  In  all  these  computations  we  assumed  that  the  exhaust  flow  is  as  in  the  typical 
HF/DF  laser  case  (Table  4-1  below),  and  that  the  ambient  density  and  velocity  are  NA=  1  x  1016 
(molecules,  m3)  and  UA  =  8  (km/ sec).  As  an  approximation  we  further  assumed  that  the  sole  ambient 
species  is  molecular  nitrogen  (molecular  weight  WA  =  28)  and  that  all  binary  collision  cross-sections 
are  uniformly  given  by  <r  =  rcD2,  where  D  is  the  molecular  diameter  (Table  4-1).  In  each  computation 
we  evaluated  the  combined  HF  +  DF  flux  by  assuming  that  the  molar  fraction  of  DF  is  zero  and  the 
molar  fraction  of  HF  is  the  combined  value  for  both  species  (Table  4-1) :  .091  +  .135  =  .226  .  This  is 
justified  by  the  relatively  small  difference  in  molecular  weight  (just  5%)  between  these  two  species. 

Three  sets  of  flux  computation  were  performed  as  follows  : 

(a)  Incidence- plane  (tp4=0)  target  points  at  various  distances  from  the  nozzle  lip  (X  =.l  to 
Xs=  10  (m)),  and  at  constant  incidence  angle  (\yA  =  20°).  The  results  are  shown  in  Fig.  4-1.  We 
observe  that  the  flux  is  fairly  insensitive  to  Xs.  Also  shown  in  Fig.  4-1  are  flux  computations 
where  the  ring-symmetric  CRW  flow  is  approximated  as  a  planar  CRW  (Prandtl- Meyer  flow), 
rather  than  the  power-law  as  in  Eq.  (2-1)  above.  The  planar  case  exhibits  a  somewhat  higher 
flux,  particularly  at  large  Xs. 

(b)  Incidence-plane  (<p4=0)  target  points  at  Xs=  I  (m)  and  at  various  incidence  angles  (y4=0  to 
\y  v  =  40°).  A  polar  representation  of  the  results  is  given  in  Fig.  4-2.  Note  the  sharp  decrease  in 
flux  as  the  incidence  angle  yA  approaches  the  plume  limiting  angle  v/f=41°. 

fc)  Azimuth  angle  variation  ((p4  =  0  to  (pA=180°)  at  a  constant  location  ( Xs  =  1  (ml)  and  at  a 
constant  angle  of  incidence  (iy4  =  20°).  A  polar  representation  of  the  results  is  given  in  Fig.  4-3. 
Observe  that  flux  becomes  sensitive  to  azimuth  angle  (p  ^  only  past  (p4  =  90°.  where  shadowing 
by  the  cylindrical  spacecraft  becomes  increasingly  dominant. 


In  addition  to  return  flux  we  also  computed  the  rms  velocity  of  the  arriving  molecules.  For  the 
target  points  in  group  (a),  the  rms  velocity  varied  between  6000  and  6600  (m  sec)  (the  higher  velocity 
at  smaller  Xs),  which  corresponds  to  a  kinetic  energy  of  about  4  (ev)  per  molecule  (HF). 
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The  maximum  return  flux  arriving  at  the  spacecraft  is  about  0.15  x  1019  (molecules/m2sec),  which 
corresponds  to  a  surface  deposition  rate  of  about  300  monolayers  (HF+DF)  per  hour.  This  level  of 
contaminating  flux  may  seem  to  be  not  outright  negligible;  however,  since  return  flux  is  proportional 
to  ambient  density,  it  will  be  scaled  down  considerably  at  higher  altitudes  (and  lower  ambient 
densities). 

We  observe  that  the  maximum  return  flux  constitutes  a  fraction  of  about  2%  of  the  incident 
ambient  flux.  This  return  flux  ratio  is  roughly  maintained  at  almost  all  target  points  and  attitude 
angles  in  groups  (a),  (b)  and  (c).  The  only  exceptions  are  incidence  angles  near  the  limiting  cone 
(i|/  =  \pf)  or  at  azimuth  angles  <pA  >  125°  where  shadowing  becomes  dominant.  This  observation  is 
interpreted  as  follows. 

Consider  the  total  solid  angle  subtended  by  the  limiting  cone  (considered  to  be  infinitely  extended 
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in  the  axial  direction)  as  viewed  from  a  target  point  (for  all  lines-of-sight  fl  pointing  outward  of  the 
cylindrical  spacecraft  surface).  It  is  independent  of  target  location  due  to  the  "self-similar"  geometry. 
During  each  flux  computation,  we  also  evaluated  the  total  solid  angle  subtended  by  that  segment  of 
the  cone  over  which  the  flux  integration  was  actually  performed  (see  Section  3.2).  It  was  found  out 
that  for  all  but  the  "shadowed"  cases  (<pA  >  125°),  this  solid  angle  constituted  a  fraction  of  86  ±  1% 
of  the  solid  angle  subtended  by  the  infinite  cone.  We  interpret  this  result  as  a  hint  that  geometrical 
"view  factors"  arising  in  the  course  of  the  flux  integration,  are  not  the  dominant  factor  in  determining 
the  2%  level  of  flux  ratio.  What  then  are  the  dominant  factors? 

For  a  possible  explanation  we  turn  to  the  flux  integration  scheme  presented  in  Section  3.2.  The 
flux  ratio  is  obtained  upon  dividing  the  integrand  in  Eq.  (3-1)  by  N  VUA  and  setting  hk=  1  (since  we 
assume  a  single  species  air).  The  major  factors  in  the  flux  ratio  integrand  appear  to  be  the  no¬ 
collision  probabilities  exp(~nik(S)|  and  exp(_nk(S)l,  and  the  post-collision  directional  distribution 
function  Pjk(S,  ~F2).  The  flux-averaged  values  of  these  functions  in  the  group  (a)  computations  were 
found  to  be  as  follows  :  P|k(S,  —  H)  =  .09  to  .10  .  r|jk(S)  =  .42  to  .54  and  nk<s)  =  .35  to  .47.  The 
flux-averaged  Mach  number  for  group  (a)  points  exhibited  a  much  larger  variation  :  between  30  and 
80.  with  the  higher  Mach  numbers  obtained  at  further  target  points. 

These  results  are  interpreted  as  follows.  The  ambient  no-collision  probability  exp|~nk(S)|  is 
sufficiently  close  to  unity,  so  that  in  an  order-of-magnitude  analysis  such  as  the  present  one.  we  may 
disregard  this  factor.  If  the  velocity  ratio  in  the  nik<S)  integral  of  Eq.  (3-1)  is  assumed  to  be  unity  (its 
average  value  for  group  (a)  points  is  about  1.4).  then  the  differential  in  the  flux  S  integration  becomes 


<rN(S)dS  =  dnik(S).  This  implies  that  the  flux  S  integration  results  in  some  average  value  of  the  only 
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remaining  factor  :  hjPjk(S,  -  C 2 ).  Since  the  Cl  integration  introduces  a  factor  of  order  unity,  the 
order-of*magnitude  estimate  for  the  arriving-to-incident  flux  ratio  is  |h;Pik(S,  —  fi)]av  •  The  value  of 
this  estimate  is  Ih.P.JS, - H)1  =  .226  x  .09 *.02,  which  is  about  equal  to  the  actual  flux  ratio  for 
target  points  in  group  (a). 

When  an  exhaust  flow  and  orbital  parameters  (velocity  and  attitude)  are  specified,  pih<s,-m 
depends  on  the  choice  of  molecular  collision  model  (we  chose  hard  spheres),  while  h-  is  uniformly 
constant.  The  foregoing  reasoning  thus  establishes  the  collision  model  as  a  significant  factor  in 
determining  ambient  scattering  flux  levels,  to  the  extent  that  Pik(S,  -  Cl)  is  sensitive  to  the  choice  of 
model. 


Table  4-1.  Typical  Operating  Conditions  of  HF/DF  Laser  Exhaust 
Mole  fractions  [H  ]  =  .091  [HF]  =  .091  [H2]  =  .104  [DF]  =  .135  [He]  =  .579 


Average  molecular  weight 

7,14 

Specific  heats  ratio 

1.54 

Stagnation  temperature  and  density 

1400  (K) 

.0075  (kg.  mJ) 

Exit  Mach  number 

4.0 

Molecular  diameter  (hard  spheres) 

2.5x10* 10  (m) 

Spacecraft  diameter 

5.0  (m) 

Nozzle  aperture 

1.0  (ml 
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5.  SPACECRAFT  CHARGING 


Spacecraft  charging  is  a  major  concern  to  spacecraft  designers,  particularly  for  missions  in  GEO 
and  to  a  lesser  extent  also  in  LEO.  The  exhaust  plume  of  an  HF'DF  laser  operating  in  the 
ionosphere  (300  to  1000  km  altitude)  may  modify  significantly  the  pre-firing  charging  pattern  of  the 
spacecraft.  Two  classes  of  effects  may  lead  to  charging  modification;  they  are: 

(a)  The  exhaust  contains  large  concentrations  of  HF  and  DF  molecules  which  are  highly 
electronegative.  They  may  be  readily  ionized  by  environmental  electrons  and  change  the  existing 
spacecraft  charging  pattern. 

(b)  When  the  spacecraft  is  oriented  obliquely  relative  to  its  orbital  velocity  and  the  ambient  plasma 
impinges  at  the  plume  boundary,  the  plume  will  cast  a  "shadow"  on  the  downstream  side, 
leading  to  a  very  dissimilar  charging  fluxes  on  the  upstream  and  downstream  halves  of  the 
spacecraft. 

The  knowledge  gained  in  analyzing  the  ambient  scattering  effect  can  be  applied  to  the  assessment 
of  the  effects  of  ionospheric  plasma  on  spacecraft  charging.  We  first  consider  the  upstream  side  of 
the  spacecraft  as  mentioned  in  (a)  above. 

We  contend  that  the  exhaust-plasma  interaction  will  not  drastically  alter  the  charging  pattern  of 

the  upstream  half.  This  assessment  is  established  as  fellows.  Consider  the  fact  that  ionospheric 

plasma  has  a  particle  number  density  no  higher  than  1012  (m*3)  and  energy  per  particle  of  at  most 

1  (ev)  (excluding  the  auroral  plasma  of  polar  zones  or  events  of  sun  storms,  where  the  energy  per 

particle  is  much  higher).  Significantlv,  the  Debve  length  at  the  highest  plasma  density  ;s  verv  small  : 

,  ’  -l 

only  about  10  (m);  the  largest  Debye  length  in  the  ionosphere  is  10  <m)  [^|  .  Ion  thermal  velocity 

is  typically  lower  than  orbital  velocity,  but  electron  velocity  is  considerably  higher  than  orbital 

velocity  (at  1  ev  the  electron  velocity  is  about  l!e  =  6x  |()~  msec).  Hence,  ions  would  typically 

impinge  at  the  plume  as  a  uniform  ion  beam  with  the  orbital  velocity  (like  ambient  molecules),  while 

electrons  are  expected  to  impinge  at  the  plume  with  their  random-oriented  thermal  velocity. 

In  view  of  the  results  of  ambient  scattering  analysis  (Ch.  3  and  4  above),  and  since  ions  are  subject 
to  similar  collision  process  with  exhaust  molecules  as  neutrals,  ions  will  be  stopped  at  the  plume 
fringes  much  like  ambient  molecules.  By  virtue  of  the  smali  Debye  length  (typically  much  smaller 
than  the  stopping  distance),  electrons  would  not  penetrate  any  further  than  ions,  regardless  of  their 


collision  cross-section  with  exhaust  molecules.  The  familiar  plasma  sheath  that  forms  at  a  solid 
surface,  is  hence  replaced  at  the  plume,  plasma  boundary-  by  a  typically  neutral  layer  whose  thickness 
is  of  the  order  of  an  ion;  neutral  mean  free  path,  but  much  larger  than  the  Debye  length.  Only  at  the 
upper  altitude  range  of  the  ionosphere  does  the  Debye  length  become  comparable  to  a  plume 
boundary  mean  free  path  (about  .1  m),  but  there  plasma  density  and  flux  are  several  orders  of 
magnitude  lower  and  charging  modification  is  not  likely  to  be  significant  at  the  relatively  short  firing 
duration  of  about  5  minutes. 

Elastically  scattered  ions  can  be  deflected  towards  the  spacecraft  as  a  result  of  elastic  collisions 
with  exhaust  molecules,  much  like  neutrals.  Referring  to  our  analysis  of  the  return-to-ambient  flux 
ratio  (Ch.  4  above),  it  is  clear  that  the  relevant  ratio  here  will  be  about  l/47t,  i.e.,  of  the  order  of  10°  o 
(this  is  due  primarily  to  the  role  played  by  the  elastic  directional  distribution  function  —  see  Ch.  4). 
A  change  in  the  plasma-to-surface  current  of  that  order  is  hence  possible,  but  unlikely  to  affect 
spacecraft  design  or  operation  significantly.  The  reason  is  that  a  design  capable  of  smoothing  away 
the  inhomogeneous  charge  flux  at  oblique  attitudes,  will  not  be  sensitive  to  a  change  in  flux  pattern 
of  the  order  of  10%  (in  other  words,  potential  differences  may  be  amplified  by  10%,  which  is  hardly 
likely  in  a  sound  design  to  bring  about  arcing  or  other  threshold  phenomena). 

Another  effect  which  may  potentially  be  significant  in  the  upstream  half  is  generation  of 
electronegative  species  (HF  ,  DF  )  by  plasma  electrons  impinging  at  the  plume.  In  the  sequel,  we 
examine  the  magnitude  of  this  effect,  concluding  that  it  is  negligible. 

•  _ 

This  estimate  is  best  done  by  considering  N  .  which  is  the  rate  of  production  of  HF  and  DF  ~ 
per  unit  volume,  at  a  typical  point  in  the  exhaust  where  local  Mach  number  is  M  =  30  'this  is 
typically  the  lowest  average  Mach  number  for  the  plume  region  where  ambient  scattering  takes  place 
~  see  Ch.  4  above).  Since  energy  is  released  by  the  electronegative  ion  formation,  the  reaction 
involves  a  third  body  as  follows  : 

HF  +e“  +  M  -»  HF"  +  M  (  M  > 

where  M  is  the  third  body  molecule.  We  assume  a  simplified  classical  kinetic  modei  for  this  reaction, 
as  follows.  The  pair  HF  M  collide  with  a  frequency  proportional  to  the  local  number  density  and  HF 
molar  fraction,  and  to  the  average  relative  velocity.  An  electronegative  ion  formation  can  occur  onh 
if  an  electron  collides  with  the  pair  durinsz  their  collision,  which  lasts  t  =D/C  .  where  C  is  the 
average  relative  pair  velocity.  Based  on  this  classical  model,  and  assuming  the  same  cross-section  lor 


electronegative  ion  formation  as  for  elastic  HF  M  collisions,  the  volume  rate  of  electronegative  tor 
generation  is  given  by  : 

N”  =  (JtD3N)  Nh(irD2UeNe)  (5-2] 

where  (ttD3N)  is  the  probability  that  a  certain  HF  or  DF  molecule  will  be  in  contact  (D  being 
molecular  hard-sphere  diameter)  with  any  other  exhaust  molecule  (whose  number  density  is  N), 
When  (ttD3N)  is  multiplied  by  hN,  where  h  is  the  HF+DF  molar  fraction  (Table  4-1),  the  combined 
term  reads  as  the  number  of  colliding  HF  M  pairs  per  unit  volume.  Assuming  the  electronegative 
formation  cross-section  is  also  trD2,  the  factor  rcD2UeNe  where  Ue  and  Ne  are  electron  velocity  and 
number  density,  renders  the  expression  for  electronegative  generation  rate  per  unit  volume.  We  note 
that  Cf  cancels  out  in  deriving  Eq.  (5-2),  so  that  N  ~  does  not  depend  on  temperature.  This  supports 
the  use  of  the  kinetic  approximation  in  regions  of  continuum  breakdown  (plume  fringes  are  such 
regions). 

•  —  •  __  ^ 

How  is  the  relative  magnitude  of  N  decided?  To  do  that  we  multiply  N  by  X  =  1/ttD‘N,  which 

is  the  mean  free  path  for  a  fast  moving  particle  that  penetrates  the  plume.  This  expression  is  justified 
by  the  fact  that  most  incident  particles  do  collide  within  a  distance  of  order  X.  and  when  the  particles 
are  plasma  ions.  electrons  will  adhere  to  ion  spatial  distribution  by  virtue  of  the  small  Debye  length 
■  smaller  than  Xi.  Thus.  XN  is  the  rate  of  electronegative  ion  generation  per  unit  area  of  plumei 
boundary.  The  ratio  p  between  this  rate  and  the  incident  electron  flux  is  : 

P"  «  XN~/N  If  =  <7TD3N)h  =  2.2  *  10'10  (5-3) 

where  N  =  2  *  IO1^  ^  m  .»  which  corresponds  to  Mach  number  ,\I  =  30  in  the  typical  case  (Table  4-1 1 
The  fraction  of  electron  flux  captured  by  HF  and  DF  exhaust  molecules  to  form  electronegative  ions 
is  so  small  i  due  to  the  pair-formation  term  <7rD3N)).  that  it  cannot  appreciably  alter  the  charging  llux 
distribution  at  the  spacecraft  surface. 

Another  possible  effect  is  the  recoil  of  HF  or  DF  that  occurs  due  to  energy  released  in  the 
electronegative  formation  reaction.  The  recoiling  species  might  conceivably  reach  the  surface  and 
contaminate  it.  The  magnitude  of  the  recoil  flux  is  certainlv  no  lareer  than  B"l  N  =1300 

-i  ~  ^  r  e 

•m  see  !.  where  we  assume  the  worst  case  ilux  :  Ne=  1 0 1 2 .  I  =6x  HF  tin  sect  which  corresponds 

to  about  1  cv  energy  per  electron.  This  llux  level  is  about  3*  10‘13  monolayers  of  HF”  and  DF~ 

per  hour,  so  that  its  contribution  to  surface  contamination  is  utterly  negligible. 
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The  second  kind  of  charging  effects  (item  (b)  above)  is  due  to  the  fact  that  the  exhaust  plume  is 
impenetrable  to  ambient  plasma  (within  a  range  of  sufficiently  small  distance  from  the  spacecraft,  so 
that  no  extensive  diluting  of  the  plume  has  taken  place).  The  downstream  half  of  the  spacecraft  in 
oblique  attitude  will  be  in  the  "shadow"  with  respect  to  incident  plasma.  As  a  first  approximation  we 
may  assume  zero  plasma  flux  at  the  shadowed  surface.  More  accurately,  this  portion  of  the 
spacecraft  will  be  subject  to  a  plasma  wake  flow.  However,  it  is  quite  difficult  to  determine  the 
charging  phenomena  that  take  place  in  such  a  wake,  as  indicated  by  a  recent  work  on  solar  sails  in 
LEO  [9]  .  Thus,  a  zero  flux  at  the  downstream  half  seems  a  practical  design  assumption. 

Can  adverse  charging  effects  occur  as  a  result  of  shadowing  the  dowmstream  half?  This  question 
can  be  discussed  only  qualitatively.  The  reason  is  that  a  quantitative  analysis  requires  a  lumped- 
circuit  model  of  the  spacecraft  external  surface  [10]  .  Since  such  a  concrete  design  is  not  available,  we 
can  only  discuss  this  question  qualitatively.  Obviously,  assuming  zero  flux  to  the  downstream  half 
during  the  envisioned  5  minutes  of  laser  firing  duration,  and  requiring  that  no  appreciable  voltages 
between  the  two  halves  will  evolve,  leads  to  the  stipulation  that  the  equivalent-circuit 
Capacitance  x  Resistance  should  be  much  smaller  than  the  firing  duration. 


6.  CONCLUDING  REMARKS 


Our  major  quantitative  conclusion  is  that  for  the  relatively  high  ambient  density  assumed 
(Na=  1  x  jo16  molecules  m3  which  represents  Sunspot  Maximum  at  about  200  km)  and  for  the 
typical  HF  DF  laser  exhaust  (Table  4-1),  the  HF+DF  flux  backscattered  by  ambient  molecules  is 
several  hundred  monolayers  per  hour.  This  flux  level  may  seem  as  not  outright  negligible.  However, 
since  ambient  scattering  flux  is  proportional  to  ambient  density,  it  will  be  scaled  down  considerably  at 
the  lower  ambient  densities  of  higher  orbital  altitudes. 

The  operational  scenario  for  HF,  DF  laser  envisions  4  or  5  minutes  total  operating  time;  hence  the 
contamination  by  ambient  scattering  may  not  be  serious  due  to  short  operating  time. 

The  effects  of  laser  exhaust  plume  on  spacecraft  charging  in  the  ionosphere  were  examined.  It  was 
concluded  that  the  rate  of  electronegative  (HF  and  DF  )  production  by  impinging  electrons  was 
negligible;  the  low  rate  is  a  consequence  of  the  assumption  that  a  third  body  is  required  to  interact 
simultaneously  with  the  HF  e  or  DF  e  pair.  No  significant  modification  of  charging  pattern  is 
anticipated.  However,  at  oblique  orbital  attitudes,  the  downstream  half  of  the  spacecraft  will  be 
shadowed  from  the  oncoming  ambient  plasma.  This  fact  has  to  be  reckoned  with  in  designing  a  ring- 
symmetric  laser  spacecraft. 


The  emphasis  in  this  work  was  on  the  method  rather  than  on  results.  The  first-collision  mode!  was 
demonstrated  to  be  simple  to  implement  in  a  code.  It  is  considerably  simpler  than  the  more  general 
and  potentially  more  accurate  Monte  Carlo  methods  commonly  used  lor  simulating  rarefied  ilows  si 
\Ve  found  out  that  the  molecular  collision  mode:  was  all  important  in  determining  the  return  ;lux 
level,  which  is  hardly  surprising  for  scattering  by  single  collision,  for  the  same  reason,  the  umi-ion 
model  would  also  be  dominant  in  a  Monte  Cario  simulation  of  the  ambient  scattering  process. 


If  and  when  a  mathematical  accuracy  of  the  first-collision  approximation  is  established  :  ir  hard- 
spheres,  it  might  be  possible  to  determine  a  reaiistic  collision  mouei  r>y  comparing  computed  resuit' 
with  measurements. 


This  accuracy  may  be  established  in  either  of  two  wavs.  One  wav  is  b\  comparison  with  accurate 
Monte  Carlo  computations  i  using  hard-spheres  collision  modei  >.  I  he  other  is  to  seek  an  estimate  u 
the  error  incurred  by  considering  just  first  collisions  and  ignoring  ail  subsequent  ones.  1  ms  might  be 
achieved  bv  accounting  for  second  collisions  in  an  extended  lirst-colIiMon  model,  presided  a  simp.dieu 
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APPENDIX  A.  DESCRIPTION  OF  AMB  CODE 


A.1  Description  of  Subroutines 

We  provide  a  list  of  the  subroutines  in  the  ambient  flux  integration  code  AMB  for  ring- symmetric 

cylidrical  spacecrafts.  Each  subroutine  is  briefly  described.  Statement  are  identified  by  the 

FORTRAN  statement  number  (columns  1  through  5). 

MAIN  PROGRAM  The  300  loop  is  intended  to  enable  several  (NCASE)  reruns  with  various  data 
in  each,  all  in  a  single  run.  Upon  calling  IN'IDATl,  parameters  depending  on  data  defined 
in  INIDAT  are  re-computed.  The  200  loop  is  over  various  XSV(N'X)  target  points.  In  the 
20  loop  the  flux  integration  begins  :  FLUXC  is  for  particle  flux  and  FLXU2C  is  for  the 
rms  of  velocity  of  return  flux  molecules.  All  the  MAX  suffixed  parameters  denote  values  at 
which  the  integrand  had  the  largest  value. 

The  actual  flux  integration  commences  at  statement  1  for  the  summation  over  strips  of 
constant  RF.  This  summation  is  terminated  when  convergence  is  attained  (to  within 
EPSR).  The  inner  loop  2  is  over  azimuth  angle  PHI.  Note  that  the  target  points  are 
generally  not  in  the  plane  of  incidence  (PHIA.XE.O).  so  that  no  symmetry  can  be  assumed 
in  the  PHI  integration,  and  it  is  performed  twice  in  order  to  cover  the  entire  range  in  PHI 
(IPAR=  1  for  PHI.GT.0.  I  PA  R  =  2  for~PHI.LT.0).  The  ilux  integration  along  the  line-ol- 
sight  is  done  by  calling  FLUX. 

INIDAT  Initialization  of  data.  There  is  no  input  file  for  this  code.  IN'IDATl  is  for  parameters 
computed  from  the  data  defined  by  calling  INIDAT. 

SOF  Stopping  routine,  called  when  an  error  is  detected.  Here  we  also  trigger  a  system  error  bv 
computing  DSQRTi-l),  in  order  to  obtain  a  calling  sequence  printout  by  the  operating 
system. 

FLUX  This  routine  calls  SU’MT  for  ilux  integration  of  one  exhaust  species  at  a  time. 

LIMIT  Here  we  compute  the  point  of  intersection  of  the  line-of-sight  with  the  leading 
characteristic  cone.  If  they  do  not  intersect,  the  distance  of  the  intersecting  point  TI.IM  is 
set  to  a  very  large  number. 
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SUMT  This  is  the  line-of-sight  integration  routine.  Runge-Kutta  scheme  is  used  (even  though  an 
explicit  integral  is  computed).  Note  that  ETAK  and  ETAIK  have  to  be  computed  through 
a  separate  integration  at  each  point  of  the  line-of-sight  integration.  The  integration  step 
DT  (T=S/RF)  is  re-adjusted  at  each  integration  step.  The  integration  is  terminated  when 
convergence  is  attained  (to  within  EPST). 

FETA  Here  the  integrand  for  the  line-of-sight  flux  integration  is  evaluated.  The  hard-spheres 
collision  model  is  used  to  determine  the  post-collision  directional  distribution  factor  PIK. 
The  flux-average  of  any  variable  (such  as  UIK**2  in  present  version),  can  be  computed  by 
summing  it  multiplied  by  flux  and  subsequently  dividing  by  the  total  arriving  flux  (see  loop 
31  in  MAIN  PROGRAM). 

PATH  IK  Here  the  molecular  thickness  ETAIK  of  the  I  exhaust  species  scattered  by  the  K  ambient 
molecule,  is  computed  by  integration  along  the  line-of-sight. 

FT  This  routine  computes  the  integrand  for  the  ETAIK  integration  in  PATH  IK. 

PATHK  The  analog  to  PATHIK  for  K  ambient  molecule.  TAU  is  the  normalized  integration 
variable  along  the  trajectory  of  the  penetrating  ambient  molecule.  Note  that 
SHADOW  =  .TRUE,  when  the  trajectory  passes  through  the  cylindrical  spacecraft  surface 
before  entering  the  fan. 

FTAU  Computes  the  integrand  for  the  ETAK  integration  in  PATHK. 

FAN  Computes  the  fan  coordinates  PS1.  XP,  YP  for  a  point  on  the  line-of-sight.  It  is  used  to 
determine  the  Mach  number  and  flow  angle  from  the  power-law  approximation  (see 
MATCH). 

FANT  Computes  the  fan  coordinates  PS  I,  XP.  YP  for  a  point  on  the  ambient  molecule  trajectory. 

HMSET  Prepares  the  vector  HMV(I)  which  is  the  value  of  the  H(M)  integral  at  a  set  of  Mach 
number  values  (equally  spaced  in  inverse  Mach  number).  This  vector  is  used  to  compute 
H(M)  for  an  arbitrary  M  (see  H1NTER),  since  this  function  is  needed  in  the  power-law 
approximation  of  flow  in  a  ring-symmetric  fan.  Subsequent  routines  MFL'NC,  HINTER. 
MATCH  and  AREAF  are  all  used  to  implement  this  approximation. 


MFUNC  Computes  the  integrand  for  the  H(M)  integration  in  HMSET. 

HINTER  Computes  H(M)  for  a  given  M,  from  HMV(I)  by  linear  interpolation.  Note  that  the 
interpolation  is  done  with  inverse  Mach  number  as  the  independent  variable. 

MATCH  Here  the  approximation  to  the  "inverse  problem"  of  finding  the  Mach  number  at  a  single 
point  in  the  ring-fan  is  implemented.  An  iteration  scheme  is  used  to  determine  the  fan 
characteristic  passing  through  the  given  point  [2]  . 

AREAF  Mach  number  is  computed  from  value  of  area  ratio  function.  Newton- Raphson  iterations 
are  used. 
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A.2  Listing  of  AMB  code 


C$0PTI0NS  LIST 

C  AMBIENT.  SCATTERING  FROM  A  RING  PLUME  BY  AMBIENT  AIR. 

IMPLICIT  REAL*8( A-H, 0-Z) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/CO , ENO, EMI , D, SIGMA, TLIM, DRO, ELO, QQ,TQ, FACT, ALOGF, 

1  DPSIO , DTMAX, DETAO , ETAL IM, XSI ,  XSF 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX, NXS, NS, NSPEC, NS1 , NS2, NTAUO, NETAO, 

1  NAMB, NCASE, I CASE, I  FAN 

COMMON  /GEOM/APF, PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 

1  CPSI1 , PSIF, SPSIF,CPSIF,TPSIF, AK, SK,CK, AO, RF,XF, YF,ZF, 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN,XS,DIST,XO,YO,ZO, 

3  DYO , DEG , PSIN , ST1 , CT 1 , OMEGX, OMEGY , OMEGZ , XSV ( 21 ) 

COMMON  /EPSIL/EPSETA,EPST,EPSR 

COMMON  /EXT REM/ TEXT (5),ETAEXT(5)» ETAKXT ( 5) , PHIEXTC  5) , 

1  PS I EXT (5),EMEXT(5), FEXT ( 5 ) , WEXT ( 5) , 

2  TMAX(5) , ETAKMXC  5) , ETAMAXC  5) , PSIMAXC5) , 

3  EMMAX(5) , FMAXC5) , 

4  RFMAX( 5) , PHI FMXC  5 ),PHIMAX(5), WMAXC5) 

COMMON  /COUNTS/ ICONTC,ICONTT, ICNTOT, ICNTMX, IQTOTC 5) , ISHAD( 5) 

COMMON  /SPEC/HAV, XC(5 ) , WC( 5) , WRCC  5) , XNAMEC  5) ,QFC(5),QDC(5), 

1  QU2C(5),FLUXC(5), OMEGAC  5) , FLXU2CC  5) , URMSCC5) 

COMMON  /AMB I EN/ ENA ,UA,PSIA,PHIA,HA(3),WA(3), 

1  UAX,UAY,UAZ, AA, BA, CA,RA,XA,YA,ZA, SHADOW 

COMMON  /POINT/XP,YP,XCOR,YCOR 
LOGICAL  SHADOW 

DIMENSION  DSUMF(5) , DSUMD(5) , DSUMAXC  5 ) , DSUMU2C5) 

NCASE=1 

DO  300  ICASE=1, NCASE 
CALL  INIDAT 

GO  TO  (301,302,303,304,305,306,307,308,309,310, 

1  311,312,313,314,315,316,317,318,319,320), I CASE 

CONTINUE 
I FAN= 1 
NXS  =  3 
XSI =0 . 1  DO 
GO  TO  399 
CONTINUE 
PHIA=20 . DO/DEG 
GO  TO  399 
CONTINUE 


301 


;o2 


;o3 


304 


305 


306 


307 


303 


309 


ilO 


311 


312 


31. 


314 


315 


GO  TO  399 
CONTINUE 
PHIA=75. DO/DEG 
GO  TO  399 
CONTINUE 
PHIA=100. DO/DEG 
GO  TO  399 
CONTINUE 
PHIA=125. DO/DEG 
GO  TO  399 
CONTINUE 
PHIA=15Q . DO/DEG 
GO  TO  399 
CONTINUE 
PHIA=175. DO/DEG 
GO  TO  399 
CONTINUE 
GO  TO  399 
CONTINUE 
GO  TO  399 
CONTINUE 
GO  TO  399 
CONTINUE 
GO  TO  399 
CONTINUE 
GO  TO  399 
CONTINUE 
GO  TO  399 
CONTINUE 
GO  TO  399 
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316  CONTINUE 
GO  TO  399 

317  CONTINUE 
GO  TO  399 

318  CONTINUE 
GO  TO  399 

319  CONTINUE 
GO  TO  399 

320  CONTINUE 
GO  TO  399 

399  CONTINUE 
PRINT  101 
101  FORMAT  ( '  1 ' ) 

C 

CALL  INDAT1 
C 

DO  200  NX=1,NXS 
XSaXSV( NX) 

C  (X0,Y0,Z0)  IS  THE  POINT  AT  WHICH  FLUX  AND  DENSITY  ARE  COMPUTED. 
C  THE  NORMAL  TO  THE  SURFACE  AT  (X0,Y0,Z0)  IS  PARALLEL  TO  Y-AXIS. 
XO  =  XS 
YO=AO 
Z0  =  0. 

DO  20  NSaNSl,NS2 
FLUXC(NS)=0. 

FLXU2C( NS) =0 . 

OMEGA(NS)=0 . 

DSUMAX(NS)=0. 

IQTOT ( NS) =0 
ISHAD(NS)=0 . 

TMAXC  NS) =-l . D  44 
ETAKMX(NS)=-1 .D  44 
PHIMAX(NS) =~1 . D  44 
PHIFMX(NS)=-1 . D  44 
WMAX( NS) =-l . D  44 
PS I MAX ( NS ) =-l . D  44 
ETAMAX(NS)=-1 .D  44 
RFMAX(NS)=-1 . D  44 
EMMAXC  NS ) =-l . D  44 
FMAX(NS)=-1 . D  44 

20  CONTINUE 
RN=RMIN 

APF=A0-0.5D0*DR0XSPSIF 
NR=  0 

1  NR=NR+1 

DR=DR0 

DR=DR0*(APF/A0) 

C  BR=DR0X(1.D0+Q.4DQ*DR0/XS)XXNR 

RF=RN+DR/2 . DO 
APF=AO+RF*SPSIF 
PHISOF=DACOS(AO/APF) 

C  DPHI0=0 . IDO 

C  NPHI=PHISOF/DPHIO+2 

DPHI=PHIS0F/D8LE(NPHI) 

DO  21  NS=NS1 ,  NS2 
DSUMF ( NS ) =0 . 

DSUMU2C  NS ) =  0 . 

DSUMDC  NS ) =0 . 

21  CONTINUE 
DOMEGR=0 . 

DO  2  NP=1 , NPHI 
DO  2  I PAR= 1 , 2 

PHIF=(DBLE(NP)-0.5D0)XDPHI 
IFCIPAR.EQ.2)  PHIF=-PHIF 
C 

CALL  FLUX 
C 

CROSS1 =QMEGY 

CR0SS2=(SPSIF)*(-0MEGX)+(-CPSIFxCPHIF)x(-0MEGY)+ 

1  (-CPSIFXSPHIF)X(-OMEGZ) 

IF(CR0SS1 . LE. 0 . ) 
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PAl 


IFCCROSS2.LE.O. ) 

1CALL  SOFC 'NORMAL  TO  LIMITING  CONE  HAS  NEGATIVE  PROJECTION  ON 
10F-SIGHT 1 ) 

D0MEGA=CR0SS2#DPHI*APF*DR/DISTXX2 

DOMEGR=DOMEGR+DOMEGA 

DO  24  NS=NS1 , NS2 

DSUMF(NS)=DSUMF(NS)+DOMEGA*QFC(NS)*CROSSl 
DSUMU2(NS)=DSUMU2(NS)+D0MEGA*QU2CCNS)*CR0SS1 
IF(DSUMAXCNS) .GT.DOMEGA*QFC(NS)*CROSSl)  GO  TO  24 
DSUMAX(NS)=DOMEGA*QFC(NS)*CROSSl 
TMAX(NS) =TEXT(NS) 

ETAKMX ( NS )=ETAKXT( NS) 

PHIMAX(NS)=PHIEXT(NS)*DEG 
PHIFMXCNS)=PHIF*DEG 
WMAX(NS) =WEXT ( NS)XDEG 
PSIMAX(NS)=PSIEXT(NS)*DEG 
ETAMAXC NS) =ETAEXT( NS ) 

RFMAX(NS)=RF 
EMMAXC  NS) =EMEXT (NS) 

FMAX(NS)=QFC(NS)*XC(NS)*QO 
»  CONTINUE 
CONTINUE 

DO  26  NS=NS1,NS2 
FLUXC(NS)=FLUXC(NS)+DSUMFCNS) 

FLXU2CC  NS) =FLXU2C<  NS)+DSUMU2(NS) 

OMEGA( NS ) =OMEGACNS)+DOMEGR 
i  CONTINUE 
RN=RN+DR 

IF(NR.LE.2)  GO  TO  1 
IF(NR.GT . 99)  GO  TO  10 
DO  27  NS=NS1 , NS2 
IFC FLUXC( NS) . EQ . 0 . )  GO  TO  27 
ERR=(DSUMF(NS)/FLUXCCNS))/DOMEGR 
IF(ERR.GT.EPSR)  GO  TO  28 
p  CONTINUE 
GO  TO  10 
!  CONTINUE 
GO  TO  1 
)  CONTINUE 

DO  31  NS=NS1,NS2 
FLUXC(NS)=XC(NS)*FLUXC(NS)XQO 

0MEGA(NS)=0MEGA(NS)/C2.D0*PAI*DC0S(PS1F/2.D0)**2) 

FLXU2C(NS)=XCCNS)*FLXU2CCNS)*Q0 

URMSCC  NS )  =0 . 

IF(FLUXCCNS) .EQ.O. )  GO  TO  31 

URMSCC  NS)=DSQRT(FLXU2C(NS)/FLUXC(NS) ) 

AVERAGE  EM  (SEE  FETA) 

URMSCC  NS ) =  FLXU2C(NS)/FLUXC(NS) 

L  CONTINUE 

PRINT  ll,NX,NR,XS,RF,DR,PHISOF*DEG 
L  FORMAT (///IX, • NX, NR , XS , RF . DR , PHI S0F= * , 21 4 , 3 D1 3 . 4 , F8 . 4 , 

1  3X, ' FLUX  AND  EXTREMA  VALUES,  ALL  SPECIES:1/) 

PRINT  12 

>  FORMAT (/IX, 1  NAME  *,*  IQTOT 1 , 1  ISHAD1, 

1  1  FMAX  1 ,’  OMEGA','  TMAX ' , 

2  '  ETAKMX','  ETAMAX','  PSIMAX', 

3  '  EMMAX','  RFMAX ' , '  PI-WMAX', 

4  '  URMSC ' , '  FLUXC  /  LOG'/) 

DO  14  NS=NS1 , NS2 

DL  F  =  0  . 

IFCFLUXCCNS) .NE.O) 

1DLF=DLOG10( FLUXC (NS) )+100 .DO+1 .D-ll 
IDLF=DLF 

DLF=DLF-DBLE(IDLF) 

PRINT  13,XNAME(NS),I0T0T(NS),ISHAD(NS), FMAX ( NS ) , OMEGA ( NS ) , 

1  TMAX(NS) , ETAKMX (NS) , ETAMAX (NS) , 

2  PSIMAX (NS) , EMMAX (NS), RFMAX (NS) , 

3  180 . DO -UMAX ( NS ) , URMSCC  NS ) , 

4  FLUXCCNS) , DLF 

5  FORMAT (1X,A6,2I6,D10.3,4F8 . 4 , 4F8 . 1 , F8 . 2 , D1 0 . 3 , ' / ' , F4 . 2 ) 

*  CONTINUE 
)0  CONTINUE 
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PRINT  102 

12  F0RMATC///1X, 'END  RING  RUN',///) 

10  CONTINUE 
STOP 
END 

SUBROUTINE  INIDAT 
IMPLICIT  REAL*8( A-H,0~Z) 

REALX8  LAMDAO , LAMDA1 
CHARACTER*8  XNAME, XNAMED 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 

1  DPSI 0 , DTMAX, DETAO,ETALIM,XSI,XSF 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, 

1  NAMB , NCASE, I CASE, I  FAN 

COMMON  /GEOM/APF, PAI ,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF, AK, SK, CK, AO , RF, XF, YF, ZF, 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMIN , RMIN,XS,DIST,XO,YO,ZO, 

3  DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , OMEGZ, XSV( 21 ) 

COMMON  /EPSIL/EPSETA,EPST,EPSR 

COMMON  /EXTREM/TEXT ( 5) , ETA EXT ( 5 ) , ETAKXT ( 5) , PHI  EXT ( 5 ) , 

1  PSI EXT ( 5 ) , EMEXT ( 5 ) , FEXT ( 5 ) , WEXT( 5 ) , 

2  TMAXC  5) , ETAKMXC  5 ) , ETAMAXC  5) , PSIMAX( 5 ) , 

3  £MMAX( 5) , FMAXC5) , 

4  RFMAX( 5 ) , PHI FMXC  5 ) , PHIMAXC  5 ) , WMAXC  5) 

COMMON  /COUNT S/ICONTC, ICONTT , ICNTOT , ICNTMX, I QTOT( 5 ) , ISHAD( 5 ) 

COMMON  /SPEC/WAV, XC(5 ) , WCC5) , WRCC5) , XNAME< 5 ) , QFCC5) , QDCC 5), 

1  QU2C( 5 ) , FLUXCC  5 ) , OMEGA ( 5 ) , FLXU2CC  5) , URMSC( 5) 

COMMON  / AMBIEN/ENA ,UA,PSIA,PHIA,HA(3),WA(3), 

1  UAX,UAY,UAZ, AA, BA, CA,RA,XA,YA,ZA, SHADOW 

COMMON  /POINT/XP,YP,XCOR, YCOR 
LOGICAL  SHADOW 

DIMENSION  XCD(5),WCD(5), XNAMEDC5) 

DATA  XCD/.091D0, .09 IDO, .104D0, .135D0, .579D0/ 

DATA  WCD/1.00DQ,20.0DO,2.00DO,21 .ODO,4,OODO/ 

DATA  XNAMED/'  H  «,'  HF  ','  H2  «,'  DF  ','  HE  '/ 

DATA  IFIRST/O/ 

IFAN=2 

PAI=4.D0*DATAN(1 .DO) 

PAI2=PAI/2.DQ 
DEG=130 . DO/PAI 
AR=8 . 3143D3 
AV=6 . 022D  26 

OMEGAC  =  0 • 5  IS  FOR  HARD  SPHERE  COLLISIONS, 

AN  AVERAGE  RECOMMENDED  VALUE  IS  ABOUT  0MEGAC=0.75 
OMEGAC-O . 5D0 
NS PEC5 5 
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NS1  =  2 
N  S  2  5  2 

DO  51  NS=1 , NSPEC 
XCC  NS ) =XCD( NS  ) 

WC(NS)=WCD(NS) 

XNAMEC  NS )= XNAMED (NS ) 

51  CONTINUE 

C  COMBINE  HF  AND  DF  MOLE  FRACTIONS  INTO  HF  FRACTION 
XC(2)=XC(2)+XC(4) 

XC ( 4 )  =  0  . 

C 

A0=2 . 5D0 
EMI =4 . DO 
RH00=0 . 0075D0 
T0=1400 . DO 
G=1 . 54D0 
D=2 . 5D-1 0 
NXS  5 1 
XSI=1 . ODO 
XSF=10 . DO 
C  AMBIENT  AIR 

ENA  =  1  .  0 0 D  16 
UA=3 . D  3 
NAMB 5 3 
WA( 1 ) =28 . DO 
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MA(2)a32 . DO 
HA(3)a16 . DO 
HACl)al .DO 
HA(2)aO. 

HAC3)a0. 

PSIAa2Q. DO/DEG 
PHIAaO . OODO/DEG 
C  INTEGRATION  PARAMETERS 
NPHI =6 
NTAUOa4 
NETAO=4 
ICNTMX=100 
RMINaO. 

DR0=0 . 10D0 
DPSI0=0 . 20D0 
DTMAX=1 . ODO 
DETA0=0 . 50D0 
ETALIM=10 . DO 
EPSTsO . 5DO 
EPSRaO . 3D0 
FACT=1.D  20 
RETURN 

CXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C  COMPUTATION  OF  DATA-DEPENDENT  PARAMETERS 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
ENTRY  INDAT1 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
AlOGFaDLOG(  FACT ) 

WAV=0. 

DO  52  NSal , NSPEC 
WAV=WAV+XCCNS)xWCCNS) 

52  CONTINUE 

DO  53  NSal, NSPEC 
WRC(NS)=WC(NS)/WAV 

53  CONTINUE 
SIGMA=PAIXDXX2 
ENO=RHOOXAV/WAV 
CO=DSQRT( GXARxTO/WAV ) 

XSV(1)=XSI 

IF(NXS.EQ.l)  GO  TO  12 

DXL  =  ( DLOG(XSF)-DLOG(XSI ) )/( DBLEC  NXS )-l . DO ) 

XLI=DLOG(XSI ) 

DO  11  NX=2 , NXS 

XSVCNX)=DEXP(XII+<DBLECNX)-1 .DO)XDXl) 

11  CONTINUE 

12  CONTINUE 
G1=(G-1 . DO )/2 . DO 
G2=(G+1.D0)/(2.D0x(G-1.D0)) 

G3=G/2.D0 

G4=(G+1.DC)/CG-1.D0) 

G5=DSQRT ( (G+l . DO )/ ( G-l . DO ) ) 

G6  =  1 . D0/(G-1 . DO ) 

G7=2.D0/(G+1 .DO) 

G8=(0 .5D0X(G+1 .D0)x*2/(G-1 .DO) )xx<l . DO/ I G+l . DO ) )x 
1  ((G+1.D0)/(G-1.D0))xx((G-1.D0)/(G+1.D0)) 

G9=CG+3.D0)/(2.D0X(G-1 .DO)) 

G1 0  =  ( 7 . DO-3 . DOXG)/( 2 . D0*( G-l . DO  ) ) 

G11=C5.D0-3.D0XG)/(2.D0X(G-1 .DO)) 

G13=(2.D0-G)/(2.DQX(G-1.D0)) 

G14=G/(2.D0X(G-1.D0)) 

G15a(G+l.D0)/(3.D0-G) 

ZETA1=G5*DATAN(DSQRT(EM1*x2-1 . D0)/G5) 

AMU1=DASIN( 1 . DO/ EMI ) 

PSI 1 =PA I 2+AMU 1 
SPSI1=DSIN(PSI1) 

CPSI1=DC0S(PSI1) 

PSIF=PAI2+AMU1+ZETA1-G5*PAI2 

SPSIFaDSIN(PSIF) 

CPSIF=DCOS(PSIF) 

TPSIF=DTAN(PSIF) 

TETA1=PSI1-AMU1 
ST1=DSIN(TETA1 ) 
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CT1*DC0SCTETA1) 

QO=ENA*UA 

LAMDA0=1 . D0/(DSQRT(2. DO)*SIGMA*ENO) 

LAMDA1=LAMDA0*(1 .DO+G1*EM1*X2)**(G6-OMEGAC+0.5DO) 
AA=DCOS(PSIA) 

BA=DSIN(PSIA)*DCOS(PHIA) 

CA=DSIN(PSIA)*DSIN(PHIA) 

UAX=-UA*AA 

UAY=-UA*BA 

UAZ=-UA*CA 

XC0R=0. 

YCOR=AO 

C 

PRINT  201 ,  NSPEC, XNAME 

201  FORMATC/1X, 'SPECIES  DATA  NSPEC=',I3/ 

1  IX, 'SPECIES  NAMES  M1(2X,  A6,2X) ) 

PRINT  202, XC 

202  FORMAT <  IX, 'MOLE  FRACTION  XC= ' , 11 C F8 . 4, 2X) ) 

PRINT  203, WC 

203  FORMAT (  IX, 'MOL  .  WEIGHT  WC= ' , 11 C F8 . 4, 2X) ) 

PRINT  21,AR,AV,WAV,G,RHOO,TO, EN0,C0,D 

21  FORMATC/1X, 'THERMODYNAMIC  DATA'/ 

1  IX, * AR, AV, WAV, GAMMA= ' , 2X, 2D14 . 5, 2F9 . 3/ 

2  IX, *RHOO,TO,ENO,CO,D=',D12.4,F8 . 0, D13 . 5, 2D12 . 4) 
PRINT  22, EMI, PSI1XDEG, PSIFXDEG, 

1  AO, LAMDAO, LAMDA1 

22  F0RMAT(/1X, 'FLOW  AND  GEOMETRY  DATA'/ 

1  IX, »EM1,PSI1,PSIF=',3F9.3/ 

2  IX, ' AO , LAMDAO , LAMDA1 = ',F9.3,2D13.4) 

PRINT  23, DPSIO, DTMAX, DETAO, ETALIM, DRO , RMIN, 

1  EPST, EPSR, 

2  NPHI , NTAUO , NETAO 

23  F0RMAT(/1X, 'INTEGRATION  DATA’/ 

1  IX, 'DPSIO, DTMAX, DETAO, ETAL IM= • , 4F9 . 4/ 

2  IX, 'DRO, RMIN, =',2013.4/ 

3  IX, *EPST,EPSR=',2D12.3/ 

4  IX, 'NPHI, NTAUO, NETAO*' ,316) 
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AMB0379 
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AMB0384 
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AMB0389 

AMB0390 

AMB0391 

AMB0392 

AMB0393 

AMB0394 

AMB0395 

AMB0396 

AMB0397 


?  *1 

■N 


PRINT  24, ENA, UA,PSIA*OEG,PHIA#DEG 
F0RMATC/1X, 'ABBREVIATED  AIR  DATA'/ 

1  IX, *ENA,UA=',2D13.4/ 

2  IX, 'PSIA,PHIA=',2F9.1) 

GO  TO  (251,252),  I  FAN 


251  CONTINUE 

PRINT  2510,  I  FAN 

2510  F0RMATC/1X, 'RING-FAN  APPROXIMATED  AS  PLANAR.  IFAN=',I4) 

GO  TO  250 

252  CONTINUE 

PRINT  2520,  IFAN 

2520  F0RMATC/1X, 'RING-FAN  APPROXIMATED  BY  MATCHED  APPROXIMATION.', 

1  4X, ' IFAN= ' , 14 ) 

250  CONTINUE 
PRINT  29 

29  F0RMAT(///1X, 'END  DATA'///) 

IFCIFIRST.EQ.0.AND.IFAN.EQ.2) 

1CALL  HMSET 

I F( I  FAN . EQ . 2 )  IFIRST=IFIRST+1 

RETURN 

END 

CSOPTIONS  LIST 

SUBROUTINE  SOF(ISTOP) 

IMPLICIT  REAL  *3 ( A-H , 0-Z ) 

CHARACTERS  ISTOP(l) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,GU , G12 , G13 , G1 4 , G1 5 . 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO , EL 0 , QO , TO , FACT , AL OGF , 

1  DPSIO, DTMAX, DETAO, ETALIM, XSI,XSF 

COMMON  /NPAR/NPHI, IPAR, NP, NR, NX, NXS, NS, NSPEC, NS1 , NS2 , NTAUO, NETAO, 

1  NAMB, NCASE, I CASE, IFAN 

COMMON  /GEOM/ APF , PAI , PAI 2 , W, SW, CW , BETA , SBETA , CBETA , PSI 1 , SPSI 1 , 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, 
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AMB0421 
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AMB  0427 
AMB0428 
AMB0429 
AMB0430 
AMB0431 
AMB0432 
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SCRIPT  A1 


3  DYO , DEG, PSIN, ST1 , CT1 , OMEGX, OMEGY, OMEGZ, X5V( 21 ) 

COMMON  /EPSIL/EPSETA, EPST, EPSR 

COMMON  /EXT REM/ TEXT ( 5) , ETA EXT ( 5) , ETAKXT ( 5) , PHI  EXT (5) , 

1  PS I EXT ( 5 ) , EMEXT ( 5) , FEXT ( 5) , NEXT ( 5) , 

2  TMAXC  5 ) , ETAKMXC 5 ) , ETAMAXC  5 ) , PSIMAXC  5 ) , 

3  EMMAXC  5) , FMAXC  5) , 

4  RFMAXC  5) , PHIFMXC  5) ,PHIMAX(5), WMAXC  5 ) 

COMMON  /SOFPR/C, DSUMF, DSUMD, T, ETA , DETA, SUM, DSUM, SUMU, DSUMU 
COMMON  /SUMS/SUMFC5) , SUMDC5) , SUMU2C5) 

COMMON  /COUNTS/ICONTC, ICONTT, ICNTOT, ICNTMX, IQTOTC  5) , ISHADC5) 

COMMON  /SPEC/WAV , XC( 5 ) , WCC 5 ) , WRCC 5 ) , XNAME( 5 ) ,  QFCC  5 ) , QDCC  5 ) , 

1  QU2CC5 ) , FLUXCC5 ) , OMEGA C  5 ) , FLXU2CC  5) , URMSCC  5 ) 

PRINT  1 , ISTOP 

FORMAT (/// 1 X, 2H**, 2X,30A4,2X,2H**,///) 

PRINT  71, NS, NP, NR, NX, ICONIC, ICONTT 

FORMAT (IX, '  NS,NP,NR,NX,IC0NTC,IC0NTT=',6I6/) 

IF(NS.GT.NSPEC)  NS  =  1 

PRINT  72, RF, PHIF*DEG, PHISOF*DEG, W»DEG, BET A* DEG 
FORMAT (/IX, 'RF,PHIF,PHIS0F,W,BETA=',D14.5,4F10.3/) 

PRINT  73,C,T,TLIM, ETA 

FORMAT (/IX, ,C,T,TLIM,ETA=',4D14.5/) 

PRINT  74,DSUM,SUM,DSUMF,SUMFCNS),SUMDCNS),QDCCNS),QFCCNS), 

1  FLUXC(NS), OMEGA C  NS ) 

FORMAT (IX, *  DSUM, SUM, DSUMF, SUMF(NS) , SUMD( NS )  =  • , 5D14 . 5/ 

1  IX, 'QDCCNS),QFCCNS),FLUXCCNS),OMEGACNS)=' ,4D14.5/) 

XX=-1 . DO 

YY=DSQRTCXX)+1  .  DO 

STOP 

END 

SUBROUTINE  FLUX 
IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /GAMA/G, G1 , G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 

1  DPSIO , DTMAX,DETAO,ETALIM,XSI,XSF 

COMMON  /NPAR/NPHI , IPAR, NP , NR, NX , NXS , NS , NSPEC, NS1 , NS2 , NTAUO , NETAO , 


COMMON 

COMMON 


1  NAMB,NCASE, ICASE, IFAN 

COMMON  /GEOM/APF, PAI,PAI2,W,SW,CW, BETA , SBETA , CBETA , PSI I , SPSI 1 , 

1  CPSIl,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,AO,RF,XF,YF,ZF, 

2  PHlSOF,PHIF,$PHIF,CPHIF, DYMIN , RMIN,XS,DIST ,XO,YO,ZO, 

3  DYO , DEG, PSIN, ST1 , CTI , OMEGX, OMEGY, OMEGZ, XSV( 21 ) 

COMMON  /EPSIL/EPSETA, EPST, EPSR 

COMMON  / EXT REM/ TEXT ( 5 ) , ETAEXTC5) , ETAKXTC 5 ) , PHI EXTC 5 ) , 

1  PSI EXTC 5) , EMEXT ( 5 ) , FEXT C  5  ) , NEXT ( 5 ) , 

2  TMAXC  5 ) , ETA KMX ( 5),ETAMAX(5),PSIMAX(5), 

3  EMMAXC  5 ) , FMAXC  5 ) , 

4  R FMAXC  5) , PHIFMXC  5 ),PHIMAXC5), UMAX C  5 ) 

COMMON  /SOFPR/C, DSUMF , DSUMD , T , ETA , DETA , SUM , DSUM , SUMU , DSUMU 
COMMON  /COUNTS/ICONTC, ICONTT, ICNTOT, ICNTMX, IQTOTC 5) , ISHADC5) 

COMMON  /SPEC/ WAV , XCC 5 ) , WCC 5 ) , WRCC 5 ) , XNAMEC 5 ) , QFCC 5 ) , QDCC 5 ) , 

1  QU2CC  5  ) , FLUXCC  5 ) , OMEGA  C  5 ) , FLXU2CC  5 ) , URMSCC  5 ) 

COMMON  /SUMS/SUMFC5) ,SUMDC5) .SUMU2C5) 

EL  0  =  SIGMA*RF*ENQ 
IFCZO.NE.O. ) 

1CALL  SOFC ’THE  SCHEME  HERE  IS  NOT  WRITTEN  FOR  ZO.NE.O.’) 
YYO=CYO-AO)/XO 
PCHECK=DATANCYYO) 

I FCPCHECK.GT.PS I F-l .D-4.0R.PCHECK.IT.-1 . D-4) 

1CALL  SOFC 'FLUX  RECEIVING  POINT  WITHIN  FAN  OR  WITHIN  SPACECRAFT') 
SPHIF=DSINCPHIF) 

CPHIF=DCOSCPHIF) 

XF’RFXCPSI F 
YF=APFXCPHIF 
ZF=APF»SPHIF 
TBETA=ZF/CYF-YO) 

BETA=DATAN(TBETA) 

IFCDABSCBETA) .GT.PAI2)  CALL  SOFC ’ BETA . GT . PAI/2 ' ) 

SBETA  =  DSINC  BETA) 

CBETA  =  DCOSC  BETA ) 

DIST  =  DSQRTC  CXF-XO )XX2+CYF-Y0 )x*2+CZF-Z0 )x*2) 

CW=  C  XF-XO )/ DI ST 
SW=DSQRTC1 ,D0-CW**2) 


IS  NOT  WRITTEN  FOR  ZO.NE.O.') 
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FILE:  AMB 


SCRIPT  A1 


W=PAI2-DATAN(CW/SH) 
OMEGXsCW 
OMEGY=SW*CBETA 
OMEGZ=SW*SBETA 
CALL  LIMIT 


DO  20  NS=NS1,NS2 
SUMF(NS)=0 . 

SUMU2C  NS ) =0 . 

SUMD(NS)=0 . 

FEXT ( NS )  =  0  . 

CALL  SUMT 

SUMF(NS)=SUM 

SUMU2(NS)=SUMU 

QFC( NS ) =$UMF( NS )/FACT 

0U2CC  NS ) sSUMU2( NS )/FACT 

FEXT(NS)=FEXT(NS)/FACT 

CALL  FAN (TEXT (NS) , PS I  EXT (NS) , PHI  EXT (NS) ) 

IF(PSIEXT(NS) .LT.PSIF-1 .D-10)  CALL  S0F( * PSIEXT( NS) . LT . PSIF’ ) 
IF(PSIEXT (NS) .GT.PSI1)  PSIEXT(NS)=PSI1 
PSIO=PSIEXT ( NS ) 


20 


T=TEXT(NS) 

CALL  MATCH( T, PSI 0 , EM, TETA ) 

EMEXT(NS)=EM 
WEXT( NS ) =W 

IQTOT(NS)=IQTOT(NS)+ICONTT 

CONTINUE 

RETURN 

END 

SUBROUTINE  LIMIT 
IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL O.QO, TO, FACT, ALOGF, 

1  DPSI 0 , DTMAX,  DETAO  , ETAL IM, XSI , XSF 

COMMON  /NPAR/NPHI , IPAR, NP, NR, NX, NXS , NS , NSPEC, NS1 , NS2, NTAUO , NETAO , 


1 


NAMB . NCASE, ICASE, IFAN 


COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW, BETA, SBETA,CBETA, PSI 1,SPSI1, 


CPSI1,P$IF.SPSIF.CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF 
PHISOF,PHIF,SPHIF,CPHIF,DYMIN, RMI N , XS , DI ST , XO , YO , ZO , 
DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , OMEGZ , XSV ( 21 ) 

COMMON  /EPSIL/EPSETA.EPST, EPSR 

COMMON  / EXT REM/ TEXT (5),ETAEXT(5)» ETAKXT ( 5 ) ,  PHI EXT(  5 ) , 

PSIEXT(5) , EMEXTC  5) , FEXT(5) , WEXT( 5 ) , 

TMAX(5) ,ETAKMX(5) ,ETAMAX(5) ,PSIMAX(5) , 

EMMAX( 5 ) , FMAX( 5 ) , 

RFMAXC  5 ) , PHI FMX( 5),PHIMAX(5), WMAX( 5 ) 


COMMON  /SPEC/WAV , XC( 5 ) . WC( 5 ) , WRC ( 5 ) , XNAME( 5 ) , Q FC ( 5 ) , QDC( 5) , 


QU2C( 5 ) , FLUXC( 5 ) , OMEGA( 5 ) , FLXU2C( 5 ) » URMSC( 5 ) 
AAA=(CW/CPSI1 )**2-l . DO 
I F( AAA . LT . 1 . D-10 )  GO  TO  1 
TPSI1=SPSI1/CPSI1 
AP1=A0+XF*TPSI1 

BBB=2 . D0*(AP1#CW#TPSI1-SW*APFX(CBETA*CPHIF+SBETA#SPHIF)) 

CCC=AP1**2-APF*#2 

DDD=BBB**2-4 . D0*AAA#CCC 

TLIM=(-BBB+DSQRT(DDD) )/(2 . DO* AAA) 

TL IM=TL IM/RF 

RETURN 

CONTINUE 

TLIM=1 . D  55 

RETURN 

END 

SUBROUTINE  SUMT 
IMPLICIT  REAL *8 ( A-H . 0-Z ) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1  G16 ,G17 , G1 8 , G1 9 , G20 

COMMON  /PAR/CO, ENO, EMI , D , S IGMA , TL IM , DRO , EL  0 , QO , T 0 , FACT , AL OGF, 

1  DPSIO.DT MAX, DETAO, ET ALIM, XSI,XSF 

COMMON  /NPAR/NPHI, I  PAR, NP, NR, NX, NXS, NS, NSPEC, NS1,NS2, NTAUO, NET AO, 
1  NAMB, NCASE, ICASE, IFAN 

COMMON  /GEOM/APF ,PAI,PAI2,W,SH,CW, BETA , SBETA , CBETA , PSI 1 , SPSI 1 , 
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SCRIPT  A1 


f 


1 

2 

3 


COMMON 

COMMON 


1 

2 

3 

4 


COMMON 

COMMON 

COMMON 


1 


CPSIl,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,AO,RF,XF,YF,ZF, 
PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, 
DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , OMEGZ, XSVC  21 ) 
/EPSIL/EPSETA,EPST,EPSR 

/EXT REM/ TEXT ( 5) , ET  A EXT ( 5 ) , ETAKXT ( 5) , PHI  EXT ( 5)  > 

PS I  EXT ( 5 ) , EMEXT ( 5 ) , FEXTC  5 ) , NEXT ( 5) , 

TMAX( 5 ) , ETAKMXC  5) , ETAMAX( 5) ,  PS I MAX ( 5 ) , 
EMMAX(5),FMAX(5), 

RFMAXC  5) ,  PHI FMX( 5 ),PHIMAX(5), WMAXC  5) 

/SOFPR/CC, DSUMF, DSUMD, T, ETA, DETA, SUM, DSUM, SUMU, DSUMU 
/COUNTS/ICONTC,ICONTT,ICNTOT,ICNTMX, IQTOTC 5 ) , ISHADC 5 ) 

/SPEC/ WAV ,  XC( 5 ) , WC( 5 ) , WRC( 5),XNAMEC5),QFC(5),QDCC5), 
QU2C(5),FLUXC(5), OMEGA<  5)»FLXU2C(5), URMSC( 5 ) 
/SUMS/SUMF C 5 ) , SUMDC 5 ) , SUMU2C  5 ) 

OF  FLUX  ARRIVING  ALONG  A  SINGLE  RAY 


10 


15 


111 


COMMON 
INTEGRATION 
DT=DPSIO 
PSIN=PSIF 
ETA1  =  0 . 

ETA3=0 . 

FETA4=0  . 

FETAU4=0 . 

T  =  0  . 

SUM  =  0  . 

SUMU=0  . 

ICONTT  =  0 

IC0NTT=IC0NTT+1 

PSIL=PSIN 

DT2=DT/2 . DO 

DT6=DT/6 . DO 

T1 =T+DT2 

T2=T+DT 

FETA1 =FETA4 

FETAU1=FETAU4 

CALL  PATHKCT1 , ETAK1 ) 

CALL  FETAC  T1 , ETA1 , ETAK1 ,GT2, FETA2, FETAU2) 

FETA3=FETA2 

CETAU3=FETAU2 

CALL  PATHKCT2, ETAK3) 

CALL  FETACT2. ETA3, ET AK3 , GT4 , FETA4 , FETAU4 ) 

DETA=DT*GT4 

D5UM=DT6*C FETA  1  +  2. DO* C FETA2+FETA3 ) +FETA4 ) 
DSUMU=DT6*CFETAUl+2. DO*C FETAU2+FETAU3 )+FETAU4 ) 

TrT+DT 
ETA=ETA3 
E” AK  =  ETAK  3 
SUM=SUM+DSUM 
'3UMU  =  3UMU  +  DSUMU 

IFCFEXTCNS) .GT.FETA4)  GO  TO  10 

FEXT ( NS )  =  FETA4 

TEXTCNS)=T 

ETAEXT C NS ) =ETA 

ETAKXT C  NS )  = ETAK 

CONTINUE 

STEP  CONTROL  ( DT  ) 

CALL  FANC  T , PSI , PHI  ) 

IFCPSI.LT. PSIF-1 . D- 1 0  )  CALL  SOFC ' PSI . LT . PSI F  '  ) 

IFCPSI.GT.PSIl)  PSI=PSI1 

P3IN=PS1 

DPS  I =P3I N-PS I L 

DTP=DT*C DPSIO/C DPSI  +  1  .  D- 1 0  ) ) 

DTE=DT*C  DETAO/C  DETA  +  1 . D-10  )  ) 

DT 1 = 1 .2D0+DT 

DT=DMIN1(DTP,DTE,DT1,DTMAX) 

IFCDT.LE.Q.)  CALL  SOFC 'COMPUTED  DT  NEGATIVE') 

CONTINUE 
IFC IPAR . LT  .  1  ) 

1  PR  I  NT  111,NR,NP,T,PSI*DEG,PHI#DEG,ETA,ETAK,SUM, DSUM/ ( SUM+1 
FORMAT C IX, ' NR, NP,T, PSI, PHI =',213,3012.3/ 

1  IX, 'ETA, ETAK, SUM, ERRR=' , 4D12. 3) 

IFC ICONTT .GT . ICNTMX) 

1CALL  SOFC 'ICONTT  TOO  LARGE’) 

IFC ICONTT . LE . 2)  GO  TO  1 


D-20) 
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FILE:  AMB 


SCRIPT  A1 


IFCETA+ETAK .GT . ETALIM)  GO  TO  100  AMB0649 

IFCT ,GT . 50 . DO  .OR.  T#RF . GT . AO )  GO  TO  100  AMB0650 

IFCSUM.EQ.O. )  GO  TO  1  AMB0651 

ERR=(DSUM/SUM)/DT  AMB0652 

IFCERR.GT.EPST)  GO  TO  1  AMB0653 

)0  CONTINUE  AMB0654 

SUM=SUM*EL0  AMB0655 

SUMU  =  SUMU#EL  0  AMB0656 

RETURN  AMB0657 

END  AMB0658 

SUBROUTINE  FETA( T, ETAIK, ETAK, GT, FET, FETU2)  AMB0659 

IMPLICIT  REAL*8(A-H,0-Z)  AMB0660 

REALMS  MU1 , MU2  AMB0661 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, AMB0662 

1  G16,G17,G18,G19,G20  AMB0663 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF,  AMB0664 

1  DPSI 0 , DTMAX, DETAO , ETAL IM, XSI , XSF  AMB0665 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB0666 

1  NAMB, NCASE, ICASE,IFAN  AMB0667 

COMMON  /GEOM/APF, PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,  AMB0668 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0669 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN,XS,DIST,X0,Y0,Z0,  AMB0670 

3  DYO,DEG,PSIN,ST1,CT1,OMEGX, OMEGY, OMEGZ, XSVC  21 )  AMB0671 

COMMON  /EPSIL/EPSETA,EPST,EPSR  AMB0672 

COMMON  /EXT REM/ TEXT (5),ETAEXT(5),ETAKXT(5), PHI  EXT ( 5) ,  AMB0673 

1  PSIEXT(5),EMEXT(5),FEXT(5), WEXT ( 5 ) ,  AMB0674 

2  TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5),  AMB0675 

3  EMMAXC  5 ) , FMAXC 5 ) ,  AMB0676 

4  RFMAXC  5),PHIFMX(5),PHIMAXC5), WMAXC  5)  AMB0677 

COMMON  /SPEC/WAV, XC( 5) , WCC 5 ) , WRCC 5 ) , XNAMEC 5 ) , QFCC 5 ) , QDCC 5 ) ,  AMB0678 

1  QU2C(5),FLUXC(5),0MEGA(5),FLXU2C(5), URMSC( 5 )  AMB0679 

COMMON  / AMBIEN/ENA ,UA,PSIA,PHIA,HA(3),WA(3),  AMB0680 

1  UAX, UAY, UAZ, AA , BA , CA , RA , XA , YA , ZA , SHADOW  AMB0681 

LOGICAL  SHADOW  AMB0682 

COMMON  /NAGESH/PIK, UIK,UIKX,UIKY,UIKZ  AMB0683 

ETAIK=0 .  AMB0684 

I F( SHADOW )  GO  TO  1  AMB0685 

K=1  AMB0636 

I=NS  AMB06S7 

CALL  FAN( T , PSI , PHI )  AMB0638 

IFCPSI.LT. PSIF-l.D-10)  CALL  SOFC ' PSI . LT . PSI F » )  AMB0689 

IFCPSI .GT.PSI1)  PSI=PSI1  AMB0690 

PSI0=PSI  AMB0691 

CALL  MATCHCT , PSIO, EM, TETA )  AMB0692 

SPSI  =  DSI N( PSI  )  AMB0693 

CPSI  =  DCOSC  PSI  )  AMB0694 

SPHI  =  DSIN( PHI  )  AMB0695 

CPHI =DCQS( PHI )  AMB0696 

ST=DSIN( TETA )  AMB0697 

CT  =  DCOSC  TETA )  AMB0698 

G0REM=1 . D0+G1*EM**2  AMB0699 

TERMN=G0REM**G6  AMB0700 

U=EM*CO/DSQRTCGOREM)  AMB0701 

UX=U*CT  AMB0702 

UY=U*ST*CPHI  AMB0703 

UZ-U*ST*$PHI  AMB0704 

COLLISION  AMB0705 

MU1=WCCI)/(WCCI)+WA(K) )  AMB0706 

MU2  =  1 . DO -MU  I  AMB0707 

UMX=MU1#UX+MU2*UAX  AMB070S 

UMY=MU1#UY+MU2#UAY  AMB0709 

UMZ=MUl*UZ+MU2#UAZ  AMB0710 

DO TUM=OMEGX*UMX+QMEGYXUMY+ OMEGZ #UMZ  AMB0711 

URX=UX-UAX  AMB0712 

URY=UY-UAY  AMB0713 

URZ=UZ-UAZ  AMB0714 

UR=DSQRT(URX*#2+URY*»2+URZ**2)  AMB0715 

DET=D0TUM#*2+(MU2#UR)##2-(UMX**2+UMY**2+UMZ##2)  AMB0716 

IFCDET.LT. 0. )  GO  TO  1  AMB0717 

DET1 =DSQRT ( DET )  AMB0718 

UIK1=-D0TUM+DET1  AMB0719 

UIK2=-D0TUM-DET1  AMB0720 


V  AMB 

V 

I 


SCRIPT  A1 


PAOl 


IFCUIK2.GT.0.)  CALI  SOFPDOUBLE  COLLISION  OPTION  NOT  PROGRAMMED 
1  YET') 

UIKSUIK1 

I FCUIK . LE . 0 . )  GO  TO  1 

UIKX=-OMEGX*UIK 

UIKY=-OMEGY*UIK 

UIKZ=-OMEGZ*UIK 

CDEL=CD0TUM+UIK)/CMU2*UR) 

IFCCDEL  . LE. 0 . )  CALL  SOFCCDEL  NEGATIVE  NOT  PROGRAMMED  YET') 
IFCCDEL-1 .D-lO.GT.l.DO) 

1CALL  SOFCCDEL  CCOSCDELTA))  CANNOT  BE  GT.l.*) 
PIK=CUIK/CMU2*UR))**2/C4.D0*PAI#CDEL) 

IF  (PIK.LT.O.)  CALL  SOF( • PIK . LT . 0 * ) 

FET=CUR/UA)*PIK/TERMN 

UREL=DSQRTCCUX-UIKX)*X2+CUY-UIKY)**2+CUZ-UIKZ)**2) 

GT=ELO*CUREL/UIK)/TERMN 

CALL  PATHIKCT, ETAIK) 

POWER=ETAIK+ETAK-ALOGF 

EFACT=0. 

IFCPOWER . LT . 60 . DO ) EFACT=DEXP(-POWER) 

FET=FET*EFACT 

FETU2=FET*UIK**2 

IFCEM.LT. 0.)  CALL  SOFC • EM . LT . 0 * ) 

FETU2=FET*EM 
RETURN 
CONTINUE 
FET  =  0 . 

FETU2=0 . 

GT  =  0 . 

RETURN 

END 

SUBROUTINE  PATHIKCTC, ETAIK) 

IMPLICIT  REAL*8(A-H,0-Z) 

REAL*8  MU1 , MU2 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
G16,G17,G18,G19,G20 


COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0,EL0,Q0, TO, FACT, ALOGF, 


DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 


COMMON  /NPAR/NPHI, IPAR, NP , NR , NX , NXS , NS , NSPEC, NS1 , NS2 , NTAUO , NETAO , 


1 


NAMB, NCASE. ICASE, IFAN 


COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA  SBETA , CBETA , PSI 1 , SPSI 1 , 


CPSI1,PSIF,SPSIF,CPSIF,TP?.F,AK,SK,CK,A0,RF,XF,YF,ZF, 
PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, 
DYO , DEG, PS  IN , ST 1 , CT1 , OMEGX , OMEGY , OMEGZ , XSV  C  21 ) 

COMMON  /EPSIL/EPSETA, EPST, EPSR 

COMMON  /EXTREM/TEXTC  5 ) , ETAEXTC  5 ) , ETAKXTC  5 ) , PHI  EXT C  5 ) , 

PS I  EXT C  5 ) , EM EXT ( 5 )  , FEXTC 5 ) , WEXT ( 5 ) , 

TMAXC  5 ) , ETAKMXC  5 ) , ETAMAXC  5),PSIMAXC5), 

EMMAXC  5 ) , FMAX ( 5 )  , 

RFMAXC  5 ) , PHI FMXC  5).PHIMAXC5), WMAXC  5 ) 


COMMON  /SPEC/ WAV , XCC5),WCC5),WRC(5),XNAMEC5),QFCC5),QDCC5), 


QU2CC5) , FLUXCC5)  ,0MEGAC5) , FLXU2CC 5 ) , URMSCC 5 ) 


COMMON  /AMBIEN/ENA, UA, PSIA, PHIA, HAC 3) , WAC 3) , 


UAX, UAY,UAZ,AA, BA, CA,RA,XA,YA,ZA, SHADOW 

LOGICAL  SHADOW 

NETA=NETAO 

DT=TC/DBLE(NETA) 

DT2=DT/2.D0 
DT6=DT/6 .DO 
GT4=0 . 

T=0 . 

ETA=0 . 

IT  =  0 
IT=IT+1 
T1=T+DT2 
T2=T+DT 
GT1 =GT  4 

CALL  FTC T1 , GT2) 

GT3=GT2 

CALL  FTC T2, GT4 ) 

DETASDT6#CGT1+2.D0#CGT2+GT3)+GT4) 

T  =  T+DT 


AMB0721 

AMB0722 

AMB0723 

AMB0724 

AMB0725 

AMB0726 

AMB0727 

AMB0728 

AMB0729 

AMB0730 

AMB0731 

AMB0732 

AMB0733 

AMB0734 

AMB0735 

AMB0736 

AMB0737 

AMB0738 

AMB0739 

AMB0740 

AMB0741 

AMB0742 

AMB0743 

AMB0744 

AMB0745 

AMB0746 

AMB0747 

AMB0748 

AMB0749 

AMB0750 

AMB0751 

AMB0752 

AMB0753 

AMB0754 

AMB0755 

AMB0756 

AMB0757 

AMB0758 

AMB0759 

AMB0760 

AMB0761 

AMB0762 

AMB0763 

AMB0764 

AM30765 

AMB0766 

AMB0767 

AMB0768 

AMB0769 

AMB0770 

AMB0771 

AMB0772 

AMB0773 

AMB0774 

AMB0775 

AMB0776 

AMB0777 

AMB0778 

AMB0779 

AMB0780 

AMB0781 

AMB0782 

AMB0783 

AMB0734 

AMB0785 

AMB0786 

AMB0787 

AMB0788 

AMB0789 

AMB0790 

AMB0791 

AMB0792 


44 


FILE:  AMB 


SCRIPT  A1 


£TA=ETA+DETA 

IF(IT.LT.NETA)  GO  TO  1 

ETAIK=ETA 

RETURN 

END 

SUBROUTINE  FT(T,GT) 

IMPLICIT  REAL*8(A-H,0-Z) 

REALK8  MU1 , MU2 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO, ELO,QO, TO, FACT, ALOGF, 

1  DPSIO , DTMAX, DETAO , ETAL IM, XSI , XSF 

COMMON  /NPAR/NPHI , IPAR, NP, NR , NX, NXS, NS, NSPEC, NS1 , NS2, NTAUO,  NETAO, 

1  NAMB,NCASE, I CASE, I  FAN 

COMMON  /GEOM/APF, PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMI N,RMIN,XS,DIST,XO,YO,ZO, 

3  DYO, DEG, PSIN, ST1 , CT1 , OMEGX, OMEGY , OMEGZ, XSVC21 ) 

COMMON  /EPSIL/EPSETA,EPST,EPSR 

COMMON  /EXTREM/TEXT ( 5 ) , ETA EXT ( 5) , ETAKXT ( 5) , PHI  EXT (5) , 

1  PSI EXT ( 5 ) , EMEXT ( 5 ) , FEXT ( 5 ) , WEXT  C  5 ) , 

2  TMAX( 5 ) , ETAKMXC  5) , ETAMAXC  5) , PSIMAXC  5) , 

3  EMMAX( 5 ) , FMAXC  5 ) , 

4  RFMAXC  5),PHIFMX(5),PHIMAX(5), WMAXC  5 ) 

COMMON  /SPEC/ WAV , XC(5),WC(5),WRC(5), XNAMEC 5) ,QFC(5),QDC(5), 

1  QU2C(5),FLUXC(5), OMEGA! 5) ,FLXU2CC5), URMSCC  5) 

COMMON  /AMBIEN/ENA, UA, PSIA,PHIA,HA(3),WA(3), 

1  UAX , UAY , UAZ, AA, BA, CA, RA, XA , YA, ZA , SHADOW 

LOGICAL  SHADOW 

COMMON  /NAGESH/PIK , UIK,UIKX,UIKY,UIKZ 
K  =  1 
I  =NS 

CALL  FAN!  T , PSI , PHI ) 

IF(PSI.LT.PSIF-l.D-lO)  CALL  SOF( 'PSI.LT.PSIF') 

IFCPSI .GT.PSI1)  PSI=PSI1 
PSIO=PSI 

CALL  MATCH(T,PSIO,EM,TETA) 

SPSI =DSIN( PSI ) 

CP SI=DCOS(PSI ) 

SPHI=DSIN(PHI) 

CPHI=DCOS( PHI ) 

ST=DSIN(TETA) 

CT=DCOS(TETA ) 

GOREM=l ,D0+G1*EM**2 
TE.RMN  =  G0REM**G6 
U=EM*CO/DSQRT  CGOREM) 

UX=U*CT 

UY=U*ST*CPHI 

'JZ-U*ST*SPHI 

UREL=DSQRT((  UX-UIKX ) #*2+ ( UY-UIKY ) X*2+( UZ-UI KZ ) *#2 ) 

GT=ELOx(UREL/UIK)/TERMN 

RETURN 

END 

SUBROUTINE  PATHK( T , ETAK ) 

IMPLICIT  REAL *3 ( A-H , Q-Z ) 

COMMON  /GAMA/G , G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,GI4,G15, 
1  GI 6 , G1 7 ,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO, EL Q,QQ, TO, FACT, ALOGF, 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 

COMMON  /NPAR/NPHI, IPAR, NP, NR, NX, NXS, NS, NSPEC, NS1,NS2, NTAUO, NETAO, 

1  NAM3, NCASE, ICASE, IFAN 

COMMON  /GEOM/APF, PAI , PAI2 , W, SW, CW, BETA , SBETA , CBETA »  PSI 1 , SPSI 1 > 

1  CP5I1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMIN , RMIN,XS,DIST,X0,Y0,Z0, 

3  DYO, DEG, PS  IN, ST  1 ,CT1, OMEGX, OMEGY, OMEGZ, XSV( 21 ) 

COMMON  /EPSIL/EPSETA, EPST, EPSR 

COMMON  / EXTREM/TEXT ( 5 ) , ETAEXT ( 5 ) , ETAKXT ( 5 ) , PHI EXT( 5 ) , 

1  PS  I  EXT ( 5 ) , EMEXT ( 5 ) , F  EXT ( 5 ) , WEXT ( 5 ) , 

2  TMAXC5), ETAKMX(S), ETAMAXC 5 ) , PSIMAXC 5 ) , 

3  EMMAXC  5 ) , FMAX! 5 ) , 

4  RFMAXI 5 ) , PHI FMX (5),PHIMAX(5), WMAX( 5 ) 

COMMON  /COUNTS/ ICONTC, ICONTT, ICNTOT, ICNTMX, I QTOT ( 5 ) , I SHADC 5 ) 


AMB0793 

AMB0794 

AMB0795 

AMB0796 

AMB0797 

AMB0798 

AMB0799 

AMB0800 

AMB0801 

AMB0802 

AMB0803 

AMB0804 

AMB0805 

AMB0806 

AMB0807 

AMB0808 

AMB0809 

AMB0810 

AMB0811 

AMB0812 

AMB0813 

AMB0814 

AMB0815 

AMB0816 

AMB0817 

AMB0818 

AMB0819 

AMB0820 

AMB0821 

AMB0822 

AMB0823 

AMB0824 

AMB0825 

AMB0826 

AMB0827 

AMB0828 

AMB0829 

AMB0830 

AMB0831 

AMB0832 

AMB0333 

AMB0834 

AMB0335 

AMB0836 

AMB0837 

AMB0838 

AMB0839 

AMB0840 

AMB0841 

AMB0842 

AMB0843 

AMB0844 

AMB0845 

AMB0846 

AMB0847 

AMB0843 

AMB0849 

AMB0850 

AMB0851 

AMB0852 

AMB0353 

AMB0854 

AMB0855 

AMB0856 

AMB0857 

AMB0858 

AMB0859 

AMB0860 

AMB0861 

AMB0862 

AMB0863 

AMB0864 


SCRIPT  A1 


COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HA(3),WA(3),  AMB0865 

1  UAX,UAY,UAZ,AA, BA, CA,RA,XA,YA,ZA, SHADOW  AMB0866 

LOGICAL  SHADOW  AMB0867 

ETAK=0.  AMB0868 

DETERMINE  POINT  OF  ENTRY  OF  AMBIENT  TRAJECTORY  TO  FAN  AMB0869 

TRF=T*RF  AMB0870 

XC=XF+TRF*OMEGX  AMB0871 

YC=YF+TRF*OMEGY  AMB0872 

ZC=ZF+TRFXOMEGZ  AMB0873 

CHECK  SHADOW  AMB0874 

SHADOW=. FALSE.  AMB0875 

EVER=BAX*2+CA**2  AMB0876 

DETS=EVERXA0XX2-(BAXZC-CAXYC)XX2  AMB0877 

IF( DETS . LE . 0 . )  GO  TO  2  AMB0878 

DET$1=DSQRT (DETS)  AMB0879 

TAU1=(-(BAXYC+CAXZC)+DETS1)/EVER  AMB0880 

IFCTAUl . GT . 0 . )  SHADOW2 . TRUE.  AMB0881 

CONTINUE  AMB0882 

IF(SHADOW)  GO  TO  10  AMB0883 

EVER1=A0+XCXTPSIF  AMB0884 

EVER2=BAXX2+CAX*2-(AAXTPSIF)X*2  AMB0885 

EVER3=BAXYC+CAXZC-AAXEVER1XTPSIF  AMB0886 

DET=EVER3XX2-EVER2X(YCXX2+ZCXX2-EVER1XX2)  AMB0887 

IFCDET.LE.O. )  AMB0888 

1CALL  SOFCNO  INTERSECTION  OF  AMB.  TRAJ .  WITH  LIMITING  CONE')  AMB0889 

DET1=DSQRT ( DET)  AMB0890 

TAUP=( -EVER3+DET1 )/EVER2  AMB0891 

TAUM=(-EVER3-DET1)/EVER2  AMB0892 

IFCTAUP . GT . 0 .  .AND.  TAUM.GT.O.)  AMB0893 

1CALL  SOFCTWO  POSITIVE  INTERSECTIONS  WITH  LIMITING  CONE. NOT  PERMITAMB0894 


1  IN  THIS  VERSION 1 ) 

TAUF=DMAX1(TAUP,TAUM) 

IFCTAUF . LE . 0 . ) 

1CALL  SOFCNO  POSITIVE  INTERSECTION  WITH  LIMITING  CONE') 
XA=XC+TAUFXAA 
YA=YC+TAUFXBA 
ZA=ZC+TAUFXCA 

RA=DSQRT(XAXX2+(DSQRT(YAXX2+ZAXX21-A0)XX2) 

TAUF=TAUF/RA 

NTAU=NTAUO 

DTAU=TAUF/DBLE(NTAU) 

ETAK=0 . 

TAU  =  0  . 

DTAU2=DTAU/2 . DO 
DTAU6=DTAU/6 .DO 
GTAU4=0 . 

ITAU=0 

ITAU=ITAU+I 

TAU1=TAU+DTAU2 

TAU2=TAU+DTAU 

GTAU1=GTAU4 

CALL  FTAU(TAU1 ,GTAU2) 

GTAU3=GTAU2 

CALL  FTAU(TAU2,GTAU4) 

DETAK=DTAU6X(GTAU1+2.D0X<GTAU2+GTAU3)+GTAU4) 

TAU=TAU+DTAU 

ETAK=ETAK+DETAK 

IF(ITAU.LT.NTAU)  GO  TO  1 

ETAK=ETAKX(SIGMAXEN0XRA) 

RETURN 

CONTINUE 

ISHADC  NS )  =  ISHADC  NS )  +  l 

RETURN 

END 

SUBROUTINE  FTAU(TAU,GTAU) 

IMPLICIT  REALX8(A-H,0-Z) 


AMB0895 

AMB0896 

AMB0897 

AMB0898 

AMB0899 

AMB0900 

AMB0901 

AMB0902 

AMB0903 

AMB0904 

AMB0905 

AM30906 

AMB0907 

AMB0908 

AMB0909 

AMB0910 

AMB0911 

AMB0912 

AMB0913 

AMB0914 

AMB0915 

AMB0916 

AMB0917 

AMB0918 

AMB0919 

AMB0920 

AMB0921 

AMB0922 

AMB0923 

AMB0924 

AMB0925 

AME0926 

AMB0927 

AMB0928 

AMB0929 

AMB0930 


COMMON  /GAMA/G, G1 ,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0931 
I  G16,G17,GI8,G19,G20  AMB0932 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF,  AMB0933 
1  DPSI 0 , DTMAX, DETAO,ETALIM,XSI,XSF  AMB0934 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB0935 
1  NAMB, NCASE, ICASE, I  FAN  AMB0936 


FILE:  AMB 


SCRIPT 


A1 


ft 

Cv 


|V. 


C 

C 

c 

c 


COMMON  /GEOM/APF, PAI , PAI2, W, SW, CW, BETA, SBETA, CBETA, PSI1 , SPSI1 , 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN,XS,DIST,X0,Y0,Z0, 

3  DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , OMEGZ, XSVC  21 ) 

COMMON  /EPSIL/EPSETA,EPST,EPSR 

COMMON  /EXTREM/TEXT ( 5) , ETAEXT C  5) , ETAKXT (5) , PHI EXTC  5) , 

1  PSIEXT(5), EMEXT C  5 ) , FEXT ( 5) , WEXT ( 5 ) , 

2  TMAX( 5 ) , ETAKMXC  5 ) , ETAMAXC  5) , PSIMAX( 5) , 

3  EMMAXC5) , FMAXC5) , 

4  RFMAX( 5 ) , PHIFMX( 5 ) , PH I MAX ( 5 ) , UMAX ( 5 ) 

COMMON  /SPEC/WAV, XC( 5 ),WC( 5) , WRCC 5) , XNAMEC 5) , QFCC 5 ) , QDCC5) , 

1  QU2CC5),FLUXCC5), OMEGAC  5) , FLXU2C( 5) , URMSCC  5) 

COMMON  /AMBIEN/ENA, UA, PSI A, PHI A, HAC  3) , WAC  3) , 

1  UAX, UAY, UAZ, AA, BA, CA, RA, XA , YA, ZA, SHADOW 

LOGICAL  SHADOW 
CALL  FANT(TAU,PSI,PHI) 

IF(PSI .LT.PSIF-1 .D-10)  CALL  SOFC ' PSI . LT . PSIF ' ) 

IFCPSI .GT.PSI1)  PSI=PSI1 
PSI0=PSI 

CALL  MATCH (T , PSIO, EM, TETA) 

SPSI=DSINC PSI ) 

CPSI=DCOSCP$I ) 

SPHI  =  DSINC  PHI ) 

CPHI=DCOSCPHI ) 

ST=DSINCTETA) 

CT=DCOSCTETA) 

GOREM=l .D0+G1XEMXK2 
TERMN=G0REM**G6 
U=EM*CO/DSQRT ( GOREM) 

UREL=DSQRTCCCT*U-UAX)*X2+CST*CPHIXU-UAY)*X2+CSTXSPHI*U-UAZ)*X2) 

GTAU=UREL/CUAXT£RMN) 

RETURN 

END 

SUBROUTINE  FAN( T , PSI , PHI) 

IMPLICIT  REAL  X8  C  A-H , O-Z) 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO , ELO , QO , TO, FACT, ALOGF, 

1  DPSIO, DTMAX.DETAO, ETALIM,XSI,XSF 

COMMON  /GEOM/APF, PAI, PAI2,W,SW,CW, BETA, SBETA, CBETA, PSI1,SPSI1, 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF, AK , SK , CK, AO , RF , XF , YF , ZF , 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMIN , RMI N , XS , DIST , XO , YO , ZO , 

3  DYO, DEG, PSIN, 3T1,CT1, OMEGX, OMEGY, OMEGZ, XSVC 21 ) 

COMMON  /POINT/XP,YP,XCOR,YCOR 

RING  FAN  GEOMETRY.  FAN  CORNER  IS  AT  C 0 , AOXCOSC PHI ) , AOXSINC PHI ) ) . 

RF  —  RADIAL  DISTANCE  ON  LIMITING  CHARACTERISTIC  OF  POINT  OF 
ENTRANCE  OF  RAY. 

DIRECTION  COSINES  OF  RAY:  OMEGX , OMEGY , OMEGZ 
TRF=T*RF 
X=XF+TRFX OMEGX 
Y=YF+TRF*OMEGY 
Z=ZF+TRF*OMEGZ 
DY=DSQRTCYXY+Z*Z)-AO 

IF(DABSCDY) . L E . 1 . D-l 0*A0 )  DY=1 . D-10*A0 
IFCDY, LT. 0. ) 

1CALL  SOFC 'POINT  X,Y,X  CANNOT  BE  CLOSER  TO  X-AXIS  THAN  RADIUS  AO') 
YY=X/DY 

PSI =PAI2-DATANCYY) 

PHI  =  DATANC  Z/Y ) 

XP=XCOR+X 

YP=AQ+DY 

RETURN 

END 

SUBROUTINE  FANTCTAU, PSI , PHI ) 

IMPLICIT  REAL #3 ( A-H , 0-Z ) 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO, EL 0,Q0, TO, FACT, ALOGF, 

1  DPSIO, DTMAX,DETAO,ETALIM,XSI.XSF 

COMMON  /GEOM/APF, PAI , PAI 2, W, 3W,CW, BETA , SBETA , CBETA , PSI 1 , SPSI 1 , 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMIN , RMI N,XS, DIST ,XO,YQ,ZO, 

3  DYO , DEG, PSIN, ST1 , CT1 , OMEGX , OMEGY , OMEGZ , XSV ( 21 ) 

COMMON  /AMBIEN/ENA, UA, PSI A, PHI  A, HAC 3 ),WA( 3), 

1  UAX, UAY, UAZ, AA, BA, CA , RA , XA , YA , ZA , SHADOW 

COMMON  /POINT/ XP,YP,XCOR,YCOR 
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SCRIPT  A1 


LOGICAL  SHADOW  AMB1009 

C  RING  FAN  GEOMETRY.  FAN  CORNER  IS  AT  ( 0 , AO*COSCPHI ) , A0*SIN(PHI ) )  .  AMB1QI0 
C  RA  --  RADIAL  DISTANCE  ON  LIMITING  CHARACTERISTIC  OF  POINT  OF  AMB1011 

C  ENTRANCE  OF  RAY.  AMB1012 

C  DIRECTION  COSINES  OF  RAY:  -AA,-BA,-CA  AMB1013 

TRA=TAU*RA  AMB1014 

X=XA-TRA*AA  AMB1015 

Y=YA-TRA*BA  AMB1016 

2=2A-TRAXCA  AMB1017 

DY=DSQRT(Y*Y+Z*Z)-AO  AMB1018 

IFCDABSCDY) . LE . 1 . D-10XA0)  DY=1.D-10*A0  AMB1019 

IFCDY.LT. 0. )  AMB1020 

1CALL  SOFC 'POINT  X,Y,X  CANNOT  BE  CLOSER  TO  X-AXIS  THAN  RADIUS  AO')  AMB1021 
YY=X/DY  AMB1022 

PSI=PAI2-DATAN(YY)  AMB1023 

PHI=DATAN(Z/Y)  AMB1024 

XP=XCOR+X  AMB1025 

YP=AO+DY  AMB1026 

RETURN  AMB1027 

END  AMBI028 

SUBROUTINE  HMSET  AMBI029 

C  SUBROUTINE  NUMBER  20  AMB1030 

IMPLICIT  REAL*8CA-H,0-Z,$)  AMB1031 

REALX8  KAPAOB, MHINV, MINVO, M,MF,M1 , M2, M3, NORM, MEXIT , LAMDOB  AMB1032 

CGMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1033 
1  G16,G17,G18,G19,G20  AMB1034 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO , ELO , QO , TO, FACT, ALOGF,  AMB1035 

1  DPSI 0 , DTMAX, DETAO,ETALIM,XSI,XSF  AMB1036 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,  AMB1037 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB1038 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO,  AMB1039 

3  DY0,DEG,PSIN,ST1,CT1,0MEGX,0MEGY,0MEGZ,XSV(21)  AMB1040 

COMMON  /GRP/ DM I NV , MHINV (101), HMV C 1 01 )  AMB1041 

COMMON  /IGRP/KHM  AMB1042 

C  A  ROUTINE  FOR  THE  C+  DERIVATIVE  DUE  TO  RING  SYMMETRY  (GRP).  AMB1043 

MEXIT=EM1  AMB1044 

KHM=51  AMB1045 

IFCKHM.GT.101)  CALL  SOFC2001  ')  AMBI046 

MINV0=1 .DO/MEXIT  AMB1047 

DMINV=MINV0/DBLE(KHM-1)  AMB1048 

M=MEXIT  AMB1049 

SUM=0 .  AMB1050 

KHMI =KHM-1  AMB1051 

DO  1  1=1, KHMI  AMB1052 

MF  =  M  AMB1053 

MHINV(I)=MINV0-D8LE(I-1)*DMINV  AMB1054 

M=1 . DO/MHINVC I )  AMB1055 

DM=M-MF  AMB105o 

Ml =M-DM  AMB1057 

M2=M-DM/2.D0  AMB1058 

M3=M  AMB1059 

CALL  MFUNC(M1,F1,ETALF1,TETA1)  AMBlObO 

CALL  MFUNC(M2,F2,ETALF2,TETA2)  AMB1061 

CALL  MFUNCCM3 , F3 , ETALF3, TETA3)  AMB1062 

SUM=SUM+DM*(Fl+4.D0*F2+F3)/6 .DO  AMB1063 

ETALF=ETALF3  AMB1064 

TETA=TETA3  AMB1365 

PSI=TETA+DASIN(1.D0/M)  AMBIO60 

NORM=( (3. D0-G)/4 . D0)*(M**2-1 . DO ) **0  . 75D0/  AMB1067 

1  (D5IN(PSI)*(1 . D0+G1*M**2)**G14)  AMB1068 

HM=SUM*NORM  AMB1069 

HMVC I ) =HM  AMB1070 

G0REM=1 . D0+G1*M**2  AMB1071 

G0R=M**2-1 . DO  AMB1072 

DELT0B  =  0  .  5DO*DSQRT(GOR)*(1 . DO/CMEXITXETAL  F)  AMB107  5 

1  +DSIN(TETA)/M)/DSIN(PSI)+G15*HM/2.D0  AMB1074 

EPSIOB=DELTOB/DSQRT(GOR)-DSIN(TETA)/(M*DSIH(PSI) )  AMB1075 

KAPAOB  =  1  .  DO  AMB1076 

IFCDABSCPAI2-TETA) .GT.l . D-6 )  AMB1077 

1KAPA0B=DTAN(TETA)#EPSI0B  AMB1078 

LAMDOB=EPSIOB-DELT03*GQREM/(GOR*DSQRT(GOR) )  AMB 107  9 

PRINT  11,I,M,HM,TETA*DEG,PSI*DEG  AMB1080 


FILE<  AMB 


SCRIPT  A1 


11  F0RMATC/1X,'  I,M,HM,TETA,PSI=*,I5,5D12.4)  AMB1081 

PRINT  12,DELTOB,EPSIOB*DEG,KAPAOB*DEG,LAMDOB*DEG  AMB1082 

12  FORMAT (  IX, ’ DELTOB, EPSI OB , KAPAOB, LAMDOB= • , 5X, 5D12 . 4 )  AMB1083 

I  CONTINUE  AMB1084 

MHINVCKHM)=0.  AMB1085 

HMV(KHM)=1 . DO  AMB1086 

RETURN  AMB1087 

END  AMB1088 

SUBROUTINE  MFUNCCM, F, ETALF, TETA )  AMB1089 

C  SUBROUTINE  NUMBER  21  AMB1090 

IMPLICIT  REAL*8(A-H,0-Z,$)  AMB1091 

REALX8  NU,  NUFUNC, M, MEXIT , MD, MDD  AMB1092 

COMMON  /GAMA/G,  Gl, G2, G3, G4, G5, G6, G7, G8,G9,G10,G11,G12,G1 3, G14,G15,AMB1 093 
1  G16,G17,G18,G19,G20  AMB1094 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF,  AMB1095 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1096 

COMMON  /GEOM/APF, PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,  AMB1097 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0» RF, XF,YF, ZF, AMB 1098 

2  PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN,XS,DIST,XO,YO,ZO,  AMB1099 

3  DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , QMEGZ, XSVC  21 )  AMB1100 

C  AMB1101 

QFCMDD) =1 . D0/DSQRT(MDD**2-1 .  DO )  AMB1102 

NUFUNCCMD) =-G5xDATAN(G5*QF(MD) l+DATANC QFCMD) )  AMB1103 

C  AMB1104 

MEXIT=EM1  AMB1105 

NU=NUFUNC(M)  AMB1106 

TETA=NUFUNC(MEXIT)+PAI2-NU  AMB1107 

G0REM=1 .D0+G1*M**2  AMB1108 

GOR=M**2-1.DO  AMB1109 

F=(MX*2)X(G0REM**G13)*DSIN(TETA)/G0RXX1 • 25D0  AMB1110 

G0REM1=1 .D0+G1*MEXIT**2  AMB1111 

G0R1=MEXIT**2-1 .DO  AMB1112 

ETALF=(CGOREM/GOREM1)X*G14)X((GOR1/GOR)X*0.25DO)  AMB1113 

RETURN  AMB1114 

END  AMB1115 

SUBROUTINE  HINTERCM, H )  AMB1116 

C  SUBROUTINE  NUMBER  22  AMB1117 

IMPLICIT  REAL*8CA-H,0-Z,$)  AMB1118 

REAL*8  MINV.M, MEXIT, MHINV  AMB1119 

COMMON  /GAMA/G , G1,G2,G3,G4,G5,G6»G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1120 
1  G16,G17,G18,G19,G20  AMB1121 

COMMON  /PAR/C0,£N0,EM1,D,SIGMA,TLIM,DR0,EL0,Q0, TO, FACT, ALOGF,  AMB1122 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1123 

COMMON  /GRP/ DMINV , MHINV ( 101 ) , HMV (101)  AMB1124 

COMMON  /IGRP/KHM  AMB1125 

C  COMPUTE  H(M)  3Y  INTERPOLATION  AMB1126 

MEXIT  =  EM1  AMB 1 127 

MINV= 1 . DO/M  AMB1128 

I=KHM-IDINT ( MI NV/ DMINV- 1 . D-9)-l  AMB1129 

IFCI.GE.l.AND.I.LT. KHM )  GO  TO  1  AMB1130 

PRINT  11, I, KHM, M, MEXIT  AMB1131 

II  FORMAT C/1X,1 I, KHM,M,MEXIT=' ,215,2014. 6/)  AMB1132 

CALL  SOFC'2201')  AMB1133 

1  CONTINUE  AMB1134 

FI  =  ( MI NV-MHI NV  C 1  +  1 ) )/ DMINV  AMB1135 

F2=l . D0-F1  AMB1136 

IFCF1.LT. -l.D-9)  CALL  SOFC'2210*)  AMB1137 

IFCF2.LT. -l.D-9)  CALL  S0F(,2211')  AMB1138 

H=F1*HMV(I)+F2*HMV(I+1 )  AMB1139 

RETURN  AMB 1 140 

END  AMB1141 

SUBROUTINE  MATCH ( T , PSI 0 , MAB , TETAAB )  AMB1142 

C  SUBROUTINE  NUMBER  23  AMB1143 

IMPLICIT  REAL*8(A-H,0-Z,$)  AMB1144 

REALMS  M, MOB, MEXIT, MAB, L AMDO B , KAPA 0 B  AMB1145 

COMMON  /GAMA/G , G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1146 
1  G16,G17,G18,G19,G20  AMB1147 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO,ELO,QO, TO, FACT, ALOGF,  AMB1148 

1  DPSIO , DTMAX, DETAO, ETALIM,XSI ,XSF  AMB1149 

COMMON  /NPAR/NPHI, I  PAR , NP , NR , NX , NXS , NS , NSPEC , NS  1 , NS2 , NTAUO , NETAO ,  AMB1150 
1  NAMB.NCASE, I CASE, I  FAN  AMB1151 

COMMON  /GEOM/APF, PAI , PAI2, W, SW, CH, BETA, SBETA,CBETA, PSI 1,SPSI1,  AMB1152 


SCRIPT  A1 


PAl 


\N 

's  AMB 


1  CPSI1 » PSIF, SPSIF, CPSIF, TPSI F, AK, SK, CK, AO, RF, XF, YF, ZF, AMB1I53 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO,  AMB1154 

3  DYO, DEG, PSIN, ST1 , CT1 ,  OMEGX, OMEGY, OMEGZ, XSVC21 )  AMB1155 

COMMON  /POINT/XP,YP,XCOR,YCOR  AMB1156 

COMMON  /GRP/DMINV,MHINV(101),HMV(101)  AMB1157 

COMMON  /IGRP/KHM  AMB1158 

MEXIT=EM1  AMB1159 

GO  TO  (101,102),IFAN  AMB1160 

101  CONTINUE  AMB1161 

C  FAN  APPROXIMATED  AS  PLANAR  AMB1162 

MAB=DSQRT ( 1 .  D0+G4/DTAN( ( PSIC  rJSIF)/G5)**2)  AMB1163 

TETAAB=PSIO-DASIN( 1 . DO/MAB)  AMB1164 

GO  TO  100  AMB1165 

102  CONTINUE  AMB1166 

C  COMPUTE  MAB  FROM  THE  INVERSE  PROBLEM  SOLUTION  AMB1167 

C0TAV=1 .DO/DTAN(PSIO)  AMB1168 

EVY=YP*DL0G(YP/YC0R)/(YP-YC0R)-1 . DO  AMB1169 

PSIN=PSIO  AMB1170 

DO  1  ITER=1,10  AMB1171 

PSI=PSIN  AMB1172 

M=DSQRT(1 .D0+G4/DTAN((PSI-PSIF)/G5)**2)  AMB1173 

M=DMAX1 (M, MEXIT )  AMB1174 

CALL  HINTERCM, HM)  AMB1175 

CALL  MFUNCCM, F, ETALF, TETA)  AMB1176 

G0REM=1 . D0+G1*MX*2  AMB1177 

G0R=M**2-1 .DO  AMB1178 

DELT0B=0.5D0*DSQRT(G0R)*(1. DO/ (MEXIT* ETALF)  AMB1179 

1  +DSIN(TETA)/M)/DSIN(PSI)+G15*HM/2.D0  AMB1180 

EPSI OB=DELTOB/DSQRT(GOR)-DSIN( TETA)/ (M*DSIN( PSI))  AMB1181 

LAMDOB=EPSIOB-DELTOB*GOREM/(GOR*DSQRT(GOR) )  AMB1182 

COTN=COTAV+LAMDOB*EVY/DSIN( PSI )  **2  AMB1183 

PSIN=PAI2-DATAN(C0TN)  AMB1184 

DPSI=PSIN-PSI  AMB1185 

IF( DABS ( DPS I ) . LT . 1 . D-6 )  GO  TO  11  AMB1186 

I  CONTINUE  AMB1187 

PRINT  12,I,ITER,PSI,PSIN,DPSI,M,XP,YP,T  AMB1188 

12  F0RMAK/1X,  ' I , ITER, PSI , PSIN, DPSI , M, XP , YP, T= '//  AMB1189 

1  IX , 214 , 7D11 . 3/ )  AMB1190 

CALL  SOFC  *  2301 ' )  AMB1191 

II  CONTINUE  AMB1192 

C  USING  MOB=M  AS  COMPUTED  FROM  THE  INVERSE  PROBLEM,  FIND  MAB.  AMB1193 

M08  =  M  AMB1194 

CALL  MFUNC(M,F, ETALF, TETA)  AMB1195 

PSI=TETA+DASIN(1 .DO/M)  AMB1196 

CALL  HINTERCM, HM)  AMB1197 

G0REM=1 . D0+G1*M**2  AMB1198 

G0R=M**2-1 . DO  AMB1199 

DELT0B=0.5D0*DSQRT(GOR)*(l. DO/(MEXIT*ETALF)  AMB1200 

1  +DSIN(TETA)/M)/DSIN(PSI )+G15*HM/2 . DO  AMB1201 

F0B=(G7*G0REM)**G2/M  AMB1202 

FAB=FOB*( YP/YCOR )**DELTOB  AMB 1203 

CALL  AREAF ( FAB , MAB  )  AMB1204 

EPSIOB=DELTOB/DSQRT(GOR)-DSIN( TETA )/(M*DSIN( PSI))  AMB1205 

KAPA0B=1 . DO  AMB1206 

IFCDABSCPAI2-TETA) .GT.l.D-8)  AMB1207 

1KAPA0B=DTAN(TETA)*EPSI0B  AMB1208 

COSTAB=DCOS( TETA )*( YP/YCOR )**(-KAPAOB)  AMB1209 

TETAAB=DACOS( COSTAB)  AMB1210 

100  CONTINUE  AMB1211 

RETURN  AMB1212 

END  AMB1213 

SUBROUTINE  AREAF( F,M)  AMB1214 

C  SUBROUTINE  NUMBER  24  AMB1215 

IMPLICIT  REAL*3(A-H,0-Z,$)  AMB1216 

REAL*8  MEXIT, MIN, M,MHINV  AMB1217 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, AMB 1218 
1  G16,G17 ,G18,G19,G20  AMB1219 

COMMON  /PAR/CO , ENO , EMI , D, SIGMA ,TLIM,DRO,ELO,QO,TO,FACT,ALOGF,  AMB1220 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1221 

COMMON  /GRP/DMINV , MHINV( 1 01 ) , HMV( 101)  AMB1222 

COMMON  /IGRP/KHM  AMB1223 

C  COMPUTE  MACH  NUMBER  M  FROM  AREA  RATIO  FUNCTION  F  AMB1224 


FILE*  AMB 


SCRIPT  A1 


C  F=(<2/(G+1))x(1+(G-1)xmxx2))x*((G+1)/(2x(G-1)>)/M 
C  INITIAL  GUESS  IS  MIN 
MEXIT=EM1 

E1=(FXMEXIT)XX(1.D0/G2)/G7 
E2=( El-1 . DO )/Gl 
E3=DMAXl(E2,MEXITxx2) 

MIN=DSQRTCE3) 

EMN=MIN 
DO  1  1=1,100 
EMO=EMN 

G0REM=1 . D0+G1XEM0XK2 
G0R=EM0xx2-l . DO 
F0=CC7XG0REM)XXG2/EM0 
DF=FO-F 

C  PRINT  123,1, EMO, EMN, FO, F, DF, GOR, GOREM 

C123  FORMAT (IX, ' I , EMO, EMN, FO, F, DF, GOR, GQREM= ' , 15, 7D12 . 4) 
DFDM=FOXGOR/(EMOXGOREM) 

DMN=DF/DFDM 
EMN=EMO-DMN 
EPSEM=DABS( DMN/EMN) 

IFCEPSEM.LT. 1 . D-l 0)  GO  TO  11 

I  CONTINUE 

CALL  SOFC ' 2401 ' ) 

II  CONTINUE 
M=EMN 
RETURN 
END 
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